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Preface 



The present Lecture Notes in Statistical Mechanics were written for the students 
of the ICTP Diploma Course in Condensed Matter Physics at the Abdus Salam 
ICTP in Trieste, Italy. 

The lectures cover classical and quantum statistical mechanics with some em- 
phasis on classical spin systems. I give also an introduction to Bose condensation 
and superfluidity but I do not discuss phenomena specific to Fermi particles, being 
covered by other lecturers. 
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Table 1 : Ensembles in statistical mechanics 



Ensemble 


Partition function 


Thermodynamic potential 


Microcanonical 


AQ(E, V, N) 


S(E, V, N) = k B In AQ(E, V, N) 


Canonical 


Z(T,V,N) 

= £• e -Ei(V,N)/k B T 

= J ASl(E, V, N)&- E{ y^l kBT dE 
Y(T,p,N) 


F(T, V, N) = -k B T In Z(T, V, N) 
G(T,p,N) = -k B T\nY(T,p,fi) 


Grand canonical 


S(T,V,|i) 

= E N ^ N/kBT Z(T,V,N) 


J(T, V, n) = F - fiN = -pV 
= -k B T\n~(T,V,fi) 
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Table 2: Thermodynamic functions 



Thermodynamic functions 
(Definition) 


Natural 
variables 


Total differential 


Entropy 

S 


(E),V,N 


dS = d(E)/T + pdV/T - fidN/T 


Internal energy 

(E) 


S,V,N 


d(E) = TdS - pdV + fidN 


Enthalpy 

H= (E) + pV 


S,p,N 


dH = TdS + Vdp + fidN 


Helmholtz free energy 
F=(E)- TS 


T,V,N 


dF = -SdT - pdV + fidN 


Gibbs free energy 

G = F + pV = Nfj, 


T,p,N 


dG = -SdT + Vdp + fidN 


Grand potential 

J = F - G = -pV 


T, V, fx 


dJ = -SdT - pdV - Ndfi 
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Table 3: Relations of different thermodynamic quantities to the partition function 



Canonical partition function 
(constant number of particles) 

Z(T, V,N) = J AQ(E, V, N)e~P E dE 

Discrete energy levels: 

Z(T,V,N)=J2 l e- m 



Free energy 

F(T, V, N) = -k B T In Z(T, V, N) G(T, P, N) = F(T, V, N) + PV 




Pressure 

P = 



dF\ 



Inverse 
compressibility 



= -V{%)t 



Entropy 

dF\ 

ar)v 



Heat capacity 
(constant V) 



a 



V 



V 




Entropy 

S = 



dG\ 

df>P 



Heat capacity 
(constant P) 



C, 



Volume 



V 



dG\ 



Isothermal 
compressibility 



y\Q P 2)T 
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Chapter 1 
Foundations 



1.1 Introduction 

A macroscopic system is, typically, composed of the order of 10 23 particles. It 
is impossible to know, investigate, or describe the exact microscopic behaviour 
of each individual particle in such a system. We must limit our knowledge to the 
average properties, thermodynamic quantities like the temperature or pressure and 
correlation functions. 

Statistical mechanics is the bridge between the microscopic and macroscopic 
world, it provides methods of calculating the macroscopic properties, like the spe- 
cific heat, from the microscopic information, like the interaction energy between 
the particles. The essential ingredient of statistical mechanics is the probability 
distribution, i.e., the collection of occupancies of different configurations (micro- 
scopic states). 

Phase Space. Consider as an example the ideal gas of N particles in a three- 
dimensional space. The Hamiltonian is 

3N 2 

^ = E|r+E% (1.1) 

i=l Zm ijjH 

(Pi is the component of the momentum, m is the particle mass, and Vy the in- 
teraction between particles. The interaction term is necessary for the system of 
particles to reach thermal equilibrium, but it can be neglected compared to the 
kinetic energy.). The 3iV coordinates and 3iV momenta form the 6iV dimensional 
phase space. Each point in the phase space represents a microscopic state of the 
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system. The microscopic state changes in time according to the canonical equa- 
tions of motion: 

9n , . OH 

Pj = - — and qj = — . (1.2) 

oqj opj 

and traces a trajectory in the phase space. 

Ensembles . The system is completely described if we know the coordinates 
and momenta of all particles (6iV variables - a terrible task). To get an average of 
a physical quantity, say A, one should make the time average over a long segment 
of the trajectory (a collection of consecutive configurations) of such a system in 
the phase space, 

1 rto+At 

(A)= lim — / A[{ qi (t), Pi (t)}]dt. (1.3) 

At^oo l\t Jt 

However, usually we know neither the exact microscopic state, i.e., the loca- 
tion of the system in the phase space, nor the trajectory. What we can know, is 
the macroscopic state (temperature, pressure, volume, ...). Instead of studying a 
particular microscopic state and the trajectory over the phase space, we investigate 
an ensemble of systems, i.e., a collection of all possible microscopic systems, all 
belonging to the same macroscopic state. 

Ergodicity. To calculate the average of a physical quantity, we substitute the 
time average over one system by an average over an ensemble of equivalent sys- 
tems at fixed time: 

(A}=[ Adq^p^pa^p^d^qd^p. (1.4) 
J w 

The integral is over the whole phase space W and p is the probability density 
in the phase space, it is the probability that a unit volume of the phase space 
is occupied. So, instead of the time average, we make an ensemble average at 
fixed time. When the time and the ensemble averages are equal, we say that the 
system is ergodic. In an alternative but equivalent definition we say that a system 
is ergodic if it evolves in such a way that it visits with equal probability all points 
of the phase space which are accessible from the initial configuration subject to 
the constraint of energy conservation. In other words, the trajectory in the phase 
space of an ergodic system will spend equal time intervals in all regions of the 
constant energy surface. Most of the physical systems of interest in statistical 
mechanics are ergodic, so we will in most of the cases replace time averages by 
the corresponding ensemble averages. 
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r(E) 



AO. 

states 



Figure 1.1: Constant-energy sur- 
face in the phase space. The area of 
the surface in the 6iV-dimensional 
phase space is T(E). 



Figure 1.2: Constant-energy shell 
in the 6iV-dimensional phase space. 
The number of states in the shell is 

AQ(E). 



For the above example of an isolated ideal gas, the total energy of the system 
is constant, H — E, the trajectory lies on the constant-energy surface T(E) in the 
phase space and the ensemble average of A is: 



(A) 



/ A({ qi:Pi })d 3N q d 3N p. 
(E) Jr(E) 



(1.5) 



T(E) Jr(E) 

(The integral is over the constant energy surface T(E).) For quantum systems as 
well as for classical systems with discrete energy levels, the number of configura- 
tions is an extremely irregular function of the energy E. Therefore we introduce 
a narrow energy interval of the width AE and the ergodic hypothesis is general- 
ized so that all the configurations, satisfying the condition that their energy is in 
the interval between E and E + AE, are equally likely to occur. The number 
of configurations Ail with the energy in the interval between E and E + AE 
is proportional to the volume of the corresponding shell in the phase space. For 
indistinguishable particles, AVL is: 

1 



AQ(E) = 



h 3N N\ 



E 



i. 



(1.6) 



E<H<E+AE 



In principle, all systems in Nature obey quantum statistical mechanics, the 
energy levels are discrete, identical particles are indistinguishable and their wave- 
functions are either symmetric (Bosons) or antisymmetric (Fermions) upon ex- 
change of two particles. Classical statistical mechanics is valid only as a special, 
limiting case when the average occupation of any single-particle quantum state is 
« 1. 
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Although we will discuss only statistical mechanics of systems in equilibrium, 
i.e., equilibrium statistical mechanics, there is a growing interest in nonequilib- 
rium statistical mechanics which investigates steady, stationary states of nonequi- 
librium systems. 



An ensemble is a collection of many configurations of a system, all describing the 
same thermodynamic state. We start with a discussion of isolated systems. 

1.2.1 Microcanonical Ensemble and the Entropy 
(Isolated Systems) 

A microcanonical ensemble represents a collection of configurations of isolated 
systems that have reached thermal equilibrium. A system is isolated from its 
environment if it exchanges neither particles nor energy with its surrounding. The 
volume, internal energy and number of particles of such a system are constant 
and are the same for all configurations that are part of the same microcanonical 
ensemble. Instead of working with constant energy surface T(E), we consider all 
the states with the total energy in a narrow energy interval between E and E+AE. 
All such configurations are lying in a thin shell in the phase space and the total 
number of configurations is Af2. 

Closely related to the ergodic hypothesis is the postulate of equal weights 
which says that in equilibrium of an isolated system, all the configurations with 
equal energy are equally probable. For a system, belonging to a microcanonical 
ensemble, therefore, the probability for being in a configuration q^Pi with the 
energy H({q i} Pi}) is: 



A£l(E) is the number of configurations having the energy in the interval E < 
7~L({ ( li,Pi}) < E + AE. All configurations {qi,Pi} which fall into this energy 
interval are occupied with equal probability p({qi,Pi}). 



1.2 Ensembles in Statistical Mechanics 




o 



E<H({qi,Pi}) <E + AE 



otherwise. 



(1.7) 
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Entropy and Probability 

Probably the most fundamental thermodynamic quantity in statistical physics is 
the entropy. Entropy is a concept, defined only for macroscopic systems. A sys- 
tem of only a few particles does not have an entropy. The entropy is a measure of 
disorder in the system, and is an additive function of state. 

Let us consider an isolated system composed of 
two subsystems, which are in contact with each other 
so that they can exchange energy and particles. As 
a result of the contact between the two subsystems, 
an additional uncertainty has been introduced into the 
system (at the beginning we knew the energy and the number of particles of each 
individual subsystem, at the end we know only the total energy and number of 
particles of the whole system). 

We postulate that the entropy is The entropy is an additive quantity, therefore: 

S = S 1 + S 2 (1.8) 

We postulate that the entropy is 




S = k B In Af2. 



(1.9) 



This relation was discovered by Ludwig Boltzmann and is therefore called the 
Boltzmann equation. The proportionality constant k B is the Boltzmann constant 
and is determined by relating the statistical definition of the entropy to the ther- 
modynamic definition. The Boltzmann equation provides the bridge between ther- 
modynamics (entropy) and statistical mechanics (number of configurations). 

The entropy plays the role of thermodynamic potential for microcanonical 
ensembles, (-S) is minimal for an isolated system in equilibrium (V, E and N 
are kept constant). From thermodynamics we know that the total differential of 
the entropy is: 

TdS(E,V,N) = dE + PdV - fidN. (1.10) 

(This is the first law of thermodynamics.) This equation provides a statistical 
definition of the temperature T, pressure P and chemical potential fx: 
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(1.12) 

N,E 



(1.13) 



E,V 



T, p and fi are intensive quantities, they are independent of N whereas the entropy 
and the energy E are extensive quantities, they are proportional to the number of 
particles N. 



1.2.2 Canonical Ensemble and the Free Energy 
(Systems at Fixed Temperature) 

A canonical ensemble is a collection of closed systems at constant temperature. 
Closed systems have fixed number of particles and volume but can exchange en- 
ergy with their surroundings. The energy of a closed system is therefore not con- 
stant. As an example we can imagine a gas of particles (air molecules, e.g.) in a 
closed container at constant temperature. 

Consider now an isolated system composed of 
a small closed (sub)system s with iV particles in 
a volume V and of a "heat bath" (b) at tempera- 
ture T. Let the two subsystems be in contact so 
that they can exchange energy but no particles. In 
equilibrium, the system s has equal temperature 
as the heat bath. Now, the energy of the small 
system, E s , is no longer fixed and the system s 
will be in any energy state E s , in each one with 
a different probability p s . The heat bath must be 
large (N b 3> N s ), so that the change in the energy 
of the bath, E b , during energy exchange between 
the small system and the bath, is negligible. Ev- 
idently, since the compound system (= system + 
bath) is isolated, the total energy E t = E b + E s is constant and the compound sys- 
tem belongs to a microcanonical ensemble for which the number of configurations 
in the energy interval between E t and E t + AE is: 




Figure 1.3: A small closed 
system (s) in contact with a 
large system, called "bath." 
s and bath can exchange en- 
ergy. 



AQt(E t ) = AQ s (E s )AQ b (E t - E s 



(1.14) 
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We first express A£l b with the entropy via the Boltzmann equation (1.9) and then 
expand In Ail b to first order in E s : 

(<99 \ E 
-^\E s = S b {E t )-^. 

(1.15) 

Thus, 

AQ b (E t - E s ) = AQ b (E t ) e"^ = C e~ pE '. (1.16) 

Here, the constant C is independent of E s and (3 is the inverse temperature, (3 = 
l/k B T. 

The (normalized) probability for the system s having the energy between E s 
and E s + dE s is: 

nfFWF - ^s(E s )An b (E t -E s )dE s 

P[ s) s ~ fdE s AQ s (E s )An b (E t -E s y {LU) 

After expressing Ail b (E t — E s ) with (1.16) we find the probability density: 

p(E s ) = ip(£ s ) e -^, (1.18) 

where p(E s ) = Ail(E s ) / AE is the density of states (number of states per unit en- 
ergy interval). Z is called the partition function and is the normalization constant, 
depending on the temperature, volume and number of particles of the system s: 



Z(T,V,N) = J p(E)e~^ E dE. (1.19) 



(Here we omit the index s.) Z(T, V, N) is a Laplace transform of the microcanoni- 
cal partition function AQ(E) = p(E)AE. For systems with discrete energy levels 
we replace the integral either by a sum over all the energy levels /: 

Z{T,V,N) = Y,9ie- pEl (1-20) 
i 

(gi is the degeneracy of the energy level /), or by a sum over all configurations i: 

Z(T,V,N) = Y, Q ~ PEi - ( L21 ) 
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The role of Z is much more important than just being the normalization con- 
stant. Z is obtained by integration over the whole phase space. In this way infor- 
mation on microscopic states of the system (coordinates and momenta of individ- 
ual particles, e.g.) becomes lost but all thermodynamic information on the system 
(with prescribed temperature, number of particles and volume) is retained and can 
be obtained from Z, as we shall see below. 

The distribution in which the particles are in equilibrium with a heat bath at 
temperature T is called a canonical distribution similarly as the ensemble defined 
by this distribution is called a canonical ensemble. For a system in equilibrium at 
a fixed temperature, the system can (at least in principle) asses any energy. The 
probability density p(E) is the same for all states with the same energy but it 
decreases exponentially with increasing energy. 



Free Energy 



In this Section we will see the relation between the average of a physical quantity 
of a canonical ensemble and the partition function Z. 

The average energy (E) and average generalized force (X) in a canonical 
distribution are given by 

(E) = ^ J E p(E)e~^ E dE = --^InZ, (1.22) 



\ I Xp(E) e^ E dE = J 



dp 

dE e -/3E dE _ 9 1 \ nZ 
dx dx p 



(1.23) 



x is the generalized coordinate, conjugate to the generalized force. Examples of 
conjugate forces and coordinates are: pressure vs. volume in gases, magnetic field 
vs. magnetization in magnetic systems; the generalized force is intensive while 
the coordinate is an extensive quantity. 

We introduce the (Helmholtz)/ree energy as: 



F(T, V, N) = -k B T In Z(T, V, N) 



(1.24) 



Like the Boltzmann equation (1.9) in case of the microcanonical ensemble, also 
this equation relates thermodynamics (free energy) and statistical mechanics (par- 
tition function). 

Relation between F and S : We take the well-known relation from thermody- 
namics: 

F(T, V, N) = (E(S, V, AO) - TS, (1.25) 
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F(T, V, N) is a Legendre transform of (E(S,V, N)), by which the independent 
variables are changed from (S, V, N) to (T, V, N). 
The total differential of F is: 

dF = -SdT - PdV + fjtdN. (1.26) 

Another Legendre transform relates the Helmholtz and the Gibbs free energies 
F(T, V, N) and G(T, P, N): 

G(T,P,N) = F(T,V,N) + PV (1.27) 

For a list of thermodynamic functions and their derivatives, see Tables 1 to 3 in 
the Appendix B. Exercise: Construct these tables for magnetic systems! 

Stability and convexity of the free energies (See Fig. 1.4). The Helmholtz free 
energy F(T, V) (for simplicity, we will omit N) is a concave (i.e, the second 
derivative is negative) function of T and convex function of V and the Gibbs free 
energy G(T, P) is always a concave function of both T and P. 




V p 



Figure 1.4: Schematic behaviour of the Helmholtz and Gibbs free energies. 

These properties follow from the fact that the system must be thermally and 
mechanically stable, therefore the specific heat and the compressibility must be 
positive quantities. We begin by noting that 

It follows that 

f) 2 G r)^ 1 
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and 

Hence, G(T,p) and F(T, V) are concave functions of T. 

Together with equation (1.26) we see that for a system in equilibrium at con- 
stant T, V, N the distribution of a system over possible configurations is such that 
F is minimal. Therefore F is also called canonical or Gibbs potential. 

The next two properties are obtained by differentiating G(T, P) and F(T, V) 
with respect to pressure and volume, respectively: 

r) 2 C r)V 
rP' F f)P 1 

(^7i) =-fe) =^7— (L32) 

Hence, G(T, P) is a concave function of pressure P and F(T, V) is a convex 
function of volume V. 

The above properties are valid also for systems that undergo phase transitions. 



1.2.3 Grand Canonical Ensemble 
(Open systems) 

Consider now a more general case with the system 
of volume V s being not only in contact with a heat 
bath at temperature T, but also exchanging parti- 
cles with the bath. As an example we can imagine 
a room at constant temperature with an open win- 
dow. In this case, neither the energy nor the num- 
ber of particles in the system (air molecules in the 
room, e.g.) are fixed. Such a system is described 
by the grand canonical distribution. 
We proceed like in case of the canonical ensemble, we consider the total sys- 
tem being part of a microcanonical ensemble, therefore the number of configura- 
tions of the total system is: 



s 

bath 



AQ(E t , N t ) = AQ^E,, N s ) AQ b (E t - E s , N t - N s ). (1.33) 
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We express Af2 b with the entropy via the Boltzmann equation and then expand 
In AQ b to first order in E s and N s : 



k B lnAQ b (E t -E s ,N t -N s ) = S b (E t - E„ N t - N s ) 

« S b (E t ,N t )-(^jE s 

= S b (E t ,N t )-f + ^ 



ON 



Thus, 



(1.34) 



AQ b (E t - E s , N t -N a ) = C e-^--^). 

Eventually we find for the probability density of a system having the energy E 
and N particles: 

p(E, N) = ^ p(E, N)q-^ e -" n \ (1.35) 

where the normalization constant S is called the grand (canonical) partition func- 
tion: 



N J 



(1.36) 



(Again, we omit the index s.) The grand partition function S is related to the 
partition function Z: 



Z(T,V,ri = Y,Z(T,V,N) e^. 



(1.37) 



N 



Thermodynamic quantities in the grand canonical ensemble 

In analogy to the free energy in the case of a canonical ensemble, we introduce 
the grand (canonical) potential J which will be used for open systems: 



J(T,V,n) = -k B T\nE 

The total differential of J is: 

dJ = -SdT - pdV - Ndp. 
The entropy can now be expressed as 



(1.38) 



(1.39) 



(1.40) 
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the pressure is: 



9V / 7> 



and the number of particles is 



\dn) Ty V d V Jt,v 

For a system in equilibrium at constant T, V, and the distribution of a system 
over possible configurations is such that the potential J is minimal. 



1.3 Analogy Between Fluids and Magnetic Systems 



Above, we have mainly discussed fluids (liquids, gases). When we want to con- 
sider, e.g., magnetic systems, the above relations between thermodynamic quanti- 
ties are easily written down for them as well. The potential energy of a magnetic 
dipole pi in a magnetic field H is 

W = -fiH (1.43) 

and the work required to increase the magnetization of a system by dM is: 

-Hdn. (1.44) 

Hence, the first law of thermodynamics for magnets reads: 

TdS = dE — HdM (1.45) 

Comparing this equation with 1.10 we see that we only have to replace the "gen- 
eralized force" - pressure P by the "force", acting on the magnetic moment, i.e., 
by the negative magnetic field: 



p -> -H. 



(1.46) 



Similarly, the "generalized coordinate" V is replaced by the "generalized coordi- 
nate" of magnetic systems, M: 



V -> M. 



(1.47) 



Notice that P and H are intensive quantities whereas V and M are extensive. 
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1.4 Fluctuations 

The microcanonical ensemble describes a collection of isolated systems where 
neither energy nor the number of particles can fluctuate. 

Energy fluctuations in the canonical ensemble: The canonical ensemble de- 
scribes a system in contact with a heat bath. For such a system, the energy is 
not exactly fixed, it can fluctuate around (E). We will now evaluate the square of 
the energy fluctuations, 

(AE) 2 = ((E - (E)) 2 } = (E 2 ) - (E) 2 . (1.48) 

In the canonical ensemble, the thermal average of the energy is 

<*> = 4| (1.49, 



and the thermal average of the energy squared is: 

l r „o .„ l d 2 z 

zW 2 

So, the square of the energy fluctuation is 



(E 2 ) = \ JE 2 P (E)e-^dE = ig. (1.50) 



{AE) ~zW~ { zW ~^ 2 ~- (L5 ) 

After expressing this with (E) (see paragraph 2.3.3) we get: 

w = -m =kBl *m (1 .52, 

(AE) 2 = k B T 2 C v . (1.53) 

Here, C v is the heat capacity, it is an extensive quantity and is proportional to N. 
The relative fluctuation of the energy is thus: 



AE k B c v T 

M = V^7' (L54) 

where e is the energy per particle and c v is the specific heat at constant volume, 
cy = Cy/N. The relative energy fluctuation is oc l/y/N and it vanishes as 
N — > oo, see Fig. 1.5. For macroscopic systems (i.e., in the limit N — > oo), the 
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Figure 1.5: Fluctuation in the energy of a small (a) and of a large (b) system. 

distribution of energies in a canonical ensemble is so sharply peaked around the 
average energy that the canonical ensemble is equivalent to the microcanonical 
ensemble. For later use we write Eq. (1.54) in the form: 



NT 



((E 2 )-(E) 2 ) (1.55) 



Fluctuations in the grand canonical ensemble: Now the number of particles is 
not fixed, so we shall calculate AiV (and AE). The calculation is analogous to 
the above calculation for the canonical ensemble. We start with the grand partition 
function 

E(T, V,n) = Y,l p{E)z-^ E ^ N) dE (1.56) 

N J 

and differentiate it twice with respect to \x by keeping j3 and V constant, 

w= mm a57) 
^=ikw- a58) 

The square of the particle number fluctuation is: 

(AN) 2 = { N* } - { Nf = {^0) tv „,9, 

For fluids at constant T and V, we obtain from the thermodynamic relation (See 
Table 1) pV = k B T InZ: 
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<9 2 ln2 V d 2 p 
dfi 2 = k^fdj^' ( ' } 

We assume that the (Helmholtz) free energy can be written as F(T, V, N) 

= Nf(T, v), where /(T, v) is the free energy per particle and depends on T and 

v = V/N. Then, 

fi(T,p) = f(T,v)+pv (1.61) 

and 



d 2 p 
dfi 2 


d /r 

dfi \v, 


)- 


1 dv 
v 2 dfi'' 


dfi 

dv 


dfi dp 
dp dv 


dp 
dv 


1 



where 



n T is the isothermal compressibility. We get finally that 



(1.62) 
(1.63) 



(N 2 ) - (Nf = {N)kBTKT . (1.64) 

v 

The square of the fluctuation in the number of particles is proportional to the 
(isothermal) compressibility! This equation shows also that AiV oc y/N and that 
AN/ (N) vanishes as N — > oo. We conclude that the grand canonical ensemble is 
equivalent to the canonical ensemble when (N) — > oo. The decision about which 
ensemble is to be used to describe the investigated system does not depend on 
whether the system is closed or open, but it is a matter of convenience. For later 
use we rewrite Eq. 1.64 in a different form: 

«r = {(N 2 ) ~ (AO 2 ) (1.65) 

The limit iV — > oo, in which the relative fluctuations become negligible, is 
called the thermodynamic limit. So far we assumed c v or k t to be constant. 
At a continuous phase transition however (see Section 1.5), these quantities and 
consequently also the corresponding fluctuations diverge on all lengthscales. This 
makes life of experimentalists and of those who make simulations on computers 
(both can deal with finite systems only) very difficult. 

The fluctuations of iV around (N) or E around (E) are - of course - time- 
dependent processes! Therefore we have to assume that the systems are ergodic 
in order to identify the time average with an ensemble average which we actually 
calculated. 
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Fluctuation-Dissipation Theorem 

Above, we have shown that (AE) 2 is proportional to the specific heat c v and that 
(AN) 2 is proportional to the compressibility. In the problem 3 we also see that 
(AM) 2 is proportional to the susceptibility. All these examples indicate that there 
is a general relation between the fluctuations of the generalized coordinate and the 
response of this coordinate to the perturbation of the conjugate generalized force. 

This is the message of the fluctuation-dissipation theorem which in a most 
general way relates the fluctuations of a physical quantity of a system in equilib- 
rium to a dissipation process which is realized when the system is subject to an 
external force driving it away from equilibrium. The theorem makes it possible 
to find the non-equilibrium properties of a system from the analysis of thermal 
fluctuations of the system in equilibrium. It is a bridge between non-equilibrium 
and equilibrium statistical mechanics. 

Now let us have a look at magnetic systems. The magnetic susceptibility, i.e., 
the response of the magnetization M to the extermal magnetic field H is: 

j ((M 2 ) - (M) 2 ) . (1.66) 
If we write M = J2i Si, the susceptibility becomes: 



X N 



X = j f T,[(^S J )-(S t )(S,)]. (1.67) 
The summand on the right-hand side is the spin-spin correlation function, 



T ij = [(S i S j )-(S i )(S j )]. (1.68) 



which describes the correlations in the spatial fluctuations of the spin around its 
thermal average. Because of the interaction between the nearest-neighbouring 
spins, one expects that the spins are correlated on short distances and that their 
fluctuations are less and less correlated as the distance increases. The correlation 
function thus decays with the distance. Often T will be assumed to have the form 

r(r) oc r- p e~ r/ «, (1.69) 

where the decay length £ is called the correlation length. As we shall see later, £ is 
a temperature-dependent quantity which plays a crucial role in phase transitions, 
it diverges at a continuous phase transition. 
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Figure 1.6: Chemical potentials of two phases in contact. When Hi = /j, 2 , the 
two phases coexist. In figure (a) the chemical potentials cross and the transition is 
first order. Dashed portions represent unstable states. In figure (b) the transition 
is continuous and the two phases meet in the critical point CP. 



1.5 Phase Transitions 



So far we have considered only homogeneous systems, i.e., single phases. In 
Nature there are, however, many systems that exist in different phases. Examples 
are solid-liquid-gas, magnetic, ferroelectric, • • •, systems. In this Section we will 
discuss transitions between such phases. 

When two phases, say 1 and 2 are in contact with each other, there is an inter- 
face separating the phases. In equilibrium, the chemical potentials of both phases 
in contact are equal, 



(1.70) 



and the two phases coexist (See Figure 1.6). 



Classification of Phase Transitions 

As one of the parameters (e.g., the temperature) is varied in such a way that the 
coexistence line is crossed, the system undergoes a phase transition between the 
two phases. At the phase transition a derivative of the free energy (or some other 
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thermodynamic potential) is discontinuous. The free energy (thermodynamic po- 
tential) becomes non-analytic in this point. The transition can be either first order 
or continuous. 

(a) First-order phase transitions. The conditions for first-order phase transi- 
tions are that the chemical potentials are equal, 



Ap = p 1 (T,p,---)-p 2 (T,p, ■■■)=(), 



(1.71) 



and that their derivatives are different, 




(1.72) 



so that the chemical potentials or the free energies cross at the transition. The 
phase with the lower chemical potential is stable and the other phase is metastable. 
An example of a first-order phase transition is the liquid-gas transition. 

(b) Continuous phase transitions. A transition is continuous when the chem- 
ical potentials of both phases are equal and when their first derivatives are also 
equal, 



A f i = fi 1 (T,p,---)- f x 2 (T,p, ■■■)=(), 




(1.73) 



(1.74) 



so that the chemical potentials or the free energies have a common tangent at the 
transition. There is no latent heat in this type of transition. 

Instead of continuous transitions, some authors write about second, third, and 
higher order transitions, depending on which derivative of the chemical poten- 
tial diverges or becomes discontinuous. This classification is somewhat arbitrary, 
therefore we will distinguish only between the continuous and first-order transi- 
tions. 



1.5.1 Order Parameter 

As an example, let us consider a simple ferromagnetic material. At high temper- 
atures, there is no magnetization and the system is rotationally invariant. At low 
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temperatures, when spontaneous magnetization occurs, the magnetization direc- 
tion defines a preferred direction in space, destroying the rotational invariance. 
We call this spontaneous symmetry breaking and it takes place at the critical tem- 
perature (for ferromagnets it is the Curie point). Since a symmetry is either absent 
or present, the two phases must be described by different functions of thermody- 
namic variables, which cannot be continued analytically across the critical point. 

Because of the reduction in symmetry, an additional parameter is needed to 
describe the thermodynamics of the low-temperature phase. This parameter is 
called the order parameter. It is usually an extensive thermodynamic variable. In 
the case of ferromagnets, the order parameter is the thermal average of the (total) 
magnetization, (M) . For the gas-liquid transition, the order parameter is either the 
volume or the density difference between the two phases. Notice that the density 
is not an extensive quantity ! Notice also that there is no obvious symmetry change 
in the gas-liquid transition. 

The basic idea of phase transitions is that near the critical point, the order 
parameter is the only relevant thermodynamic quantity. 

Some examples of the order parameters and conjugate fields are given in the 
following table: 



Systems 


Order Parameter 


Conjugate Field 


Gas-liquid 


Vg-Vl 


P-Pc 


Ferromagnets 


M 


H 


Antiferromagnets 


-^■staggered 


H staggered 


Superfluids 


(Condensate wave fn.) 


Not physical 


Superconductors 


A 

(Gap parameter) 


Not physical 



1.5.2 Critical Point 

1. Liquid-gas transitions. In Figure 1.7(a) a typical phase diagram of a solid- 
liquid-gas system in the P—T plane at constant V and N is shown. Here we would 
like to concentrate on the liquid-gas transition. The liquid and the gas phases 
coexist along a line where the transition between these two phases is first order. 
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If T is increased, the coexistence line terminates in the critical point, where the 
transition between the two phases is continuous. Beyond the critical point, there 
is no transition between the two phases, there is no singularity in any physical 
quantity upon going from one phase to the other. It is also instructive to consider 
the isotherms in the P — p plane at constant T and N, see Fig. 1.7(b). At low 
temperature, there is a large difference between the gas and liquid densities, pa 
and pl, but as the critical temperature is approached, this difference tends to zero. 
The existence of a quantity that is non-zero below the critical temperature and zero 
above it, is a common feature of critical points in a variety of different physical 
systems. p L — pais the order parameter for the liquid-gas critical point. 

2. Ferromagnets. A typical phase diagram of a ferromagnet is shown in Fig. 
1.8(a). The order parameter is the magnetization M, i.e., the total magnetic mo- 
ment. At low temperature, M is ordered either "up" or "down", depending on the 
external field H. When H changes sign, the system undergoes a first-order phase 
transition (bold line) in which the magnetization changes its direction discontin- 
uously [Fig. 1.8(b)]. As T is increased, the critical point is reached, at which 
M changes continuously upon variation of H but its derivative with respect to H, 
i.e., the susceptibility diverges. At still higher T, there is no phase transition, the 
spins are simply disordered, this is the paramagnetic region. 



1.5.3 Critical Exponents 

In this Section the behaviour of the thermodynamic quantities near the critical 
point is discussed. We suppose that very near the critical point, any thermody- 
namic quantity can be decomposed into a "regular" part, which remains finite (not 
necessarily continuous) and a "singular part" that can be divergent or have diver- 
gent derivatives. The singular part is assumed to be proportional to some power 
of \T-T C \ (orp- Pc ). 

Definition of the critical exponents. For clarity, we will define the critical ex- 
ponents on a ferromagnet. However, their definition will be general and will be 
valid for arbitrary system that undergoes a continuous phase transition. The most 
important critical exponents are: 

P = Order parameter exponent. When H = 0, the magnetization M in a 
ferromagnet is a decreasing function of T and vanishes at T c (See Fig. 1.9). 

For T very close to (but below) T c , M can be written as 

M~(T C -Tf. (1.75) 



21 



I. Vilfan 



Statistical Mechanics 



Chapter 1 




Figure 1.9: Order parameter vs. temperature and field. 

(3 is the order parameter exponent which describes the behaviour of M near Tq. 
Here and later, ~ means "the singular part is proportional to". 

5 = Exponent of the critical isotherm. When T = T c , M vanishes for H = 
but increases very rapidly if H increases. For small H, the singular part of M is 
proportional to: 

M ~ H 1/s , (1.76) 

where 5 > 1. 

7 = Susceptibility exponent. As T c is approached either from above or from 
below, the susceptibility diverges. The divergence is described with the exponent 

7: 

X ~\T-T c \-\ (1.77) 

Notice that the exponent 7 is the same for T > T c and T < T c whereas the 
proportionality constants are different. 

a = Specific heat exponent. The specific heat at constant volume (for ferro- 
magnetic systems one has to take H = 0) has a singularity at Tq. This singularity 
is described by the exponent a: 

C V ~\T-T c \~ a . (1.78) 

As usually, the exponent a is the same for T > T c and T < T c whereas the 
proportionality constants are different. 

v = Correlation length exponent. As already mentioned, the correlation length 
£ is a temperature- dependent quantity. As we shall see later, £ diverges at the 
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critical point with the exponent v. 

i~\T-T c \~ v . (1.79) 

i] = Correlation function exponent at the critical isotherm. At Tq, £ diverges 
and the correlation function T decays as a power law of r, 

T ~ r d ' 2 +\ (1.80) 

where d is the dimensionality of the system. 

Scaling. As the name suggests, the scaling has something to do with the 
change of various quantities under a change of length scale. Important scaling 
laws concerning the thermodynamic functions can be derived from the simple 
(but strong) assumption that, near the critical point, the correlation length £ is the 
only characteristic length of the system, in terms of which all other lengths must 
be measured. This is the "scaling hypothesis." 

With the aid of the scaling hypothesis we are able to get relations between the 
critical exponents. The most common scaling relations are: 



7 = P(5 - 


1) 


7 = u{2 - 


v) 


a + 2(3 + 7 


= 2 


vd = 2 - 


a 



(1.81) 



so that only two of the above exponents are independent. The last four equations 
are called the scaling laws. We shall prove some of the scaling laws later. 

Universality. Some quantities depend on the microscopic details of the system. 
Examples of such quantities are: 

- The critical temperature depends on the strength of interaction. 

- The proportionality constants (e.g., for the temperature dependence of the 
magnetization) also depend on the microscopic properties of the system. 

On the other hand, the critical exponents and the functional form of the equa- 
tion of state, of the correlation functions, etc. are independent of many details of 
the system. They depend only on global features, such as: 

- The spatial dimensionality, 

- symmetry (isotropic, uniaxial, planar,...), ... . 

Systems that have the same critical exponents and the same functional form 
of the equation of state, etc, belong to the same universality class. For example, 
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the Ising model in d = 2 forms one universality class, whereas in d = 3 it forms 
another universality class, the Ising models in d = 2 and d = 3 have different 
critical exponents because they have different spatial dimensionality. However, on 
a first glance very different physical systems often belong to the same universality 
class. Example: binary alloys and some magnetic systems can be described with 
mathematically similar models and they belong to the same universality class. 
Roughening of some crystal surfaces at high temperature (but below the melting 
temperature) and the d = 2 Coulomb gas are members of another universality 
class. 
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Chapter 2 
Classical Models 



In this Chapter we will discuss models that can be described by classical statistical 
mechanics. We will concentrate on the classical spin models which are used not 
only to study magnetism but are valid also for other systems like binary alloys or 
lattice gases. 

In principle, all particles obey quantum statistical mechanics. Only when the 
temperature is so high and the density is so low that the average separation be- 
tween the particles is much larger than the thermal wavelength A, 

/V\ 1/3 ( h 2 \ 1/2 

(-) » A, where A =(^^j , (2.1) 

the quantum statistical mechanics reduces to classical statistical mechanics and we 
can use a classical, Maxwell-Boltzmann distribution function (instead of Fermi- 
Dirac or Bose-Einstein). 

For magnetic materials, the situation is somewhat different. The Hamiltonian 
of a magnetic system is a function of spin operators which can usually not be 
directly approximated by classical vectors. Quantum models in d dimensions can 
be mapped onto classical models in d! — (d + 1) dimensions. The classical limit 
of the Heisenberg model, however, can be constructed for large eigenvalues of 
the spin operator by replacing the spin operators by three-dimensional classical 
vectors. 

There is another quantum-mechanical effect we must discuss. In quantum me- 
chanics identical particles are indistinguishable from each other whereas in classi- 
cal mechanics they are distinguishable. Therefore, the partition function of a sys- 
tem of non-interacting, non-localized particles (ideal gas!) is not just the product 
of single-particle partition functions, as one would expect from classical statistical 
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mechanics, but we must take into account that the particles are indistinguishable. 
We must divide the total partition function by the number of permutations between 
iV identical particles, AM, 

Z N = j^(Z 1 ) N , (2.2) 

where Z 1 is the one-particle partition function. The factor AM is a purely quantum 
effect and could not be obtained from classical statistical mechanics. 

In (insulating) magnetic systems we deal with localized particles (spins) and 
permutations among particles are not possible. Therefore we must not divide the 
partition function by (AM). 

A Hamiltonian of a simple magnetic system is: 

1 ^ — > — * ^ — ^ ^* 

ft — — ^ / , JjjSjSj ~ H y ^ Sj. (2.3) 

The first term is the interaction energy between the spins, J it j being the (ex- 
change) interaction energy between the spins at the sites i and j. The sum is 

over all bonds between the spins. S is the spin operator; in general it is a vec- 
tor. The order parameter is equal to the thermal and quantum-mechanical av- 
erage of S, rn = ((S)). The thermal average of the first term in Eq. (2.3) is 
E(S, m) = (E(S, {S})}. (Here, S is the entropy and {S} is a spin configuration. 
The thermal average has to be taken over all spin configurations.) E(S, m) is the 
internal energy of a magnetic system. 

The second term in (2.3) is the interaction with external magnetic field H. The 
thermal average of this term corresponds to PV in fluids. We see that (H(S, H)) 
is a Legendre transformation of E(S, m) and is thus equivalent to the enthalpy of 
fluids. 



2.1 Real Gases 

In the ideal gas we neglected the collisions between the particles. In reality, the 
atoms (molecules) have a finite diameter and when they come close together, they 
interact. The Hamiltonian of interacting particles in a box is: 



-Q 

2m 



i>j 



(2.4) 
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u{r-i — Tj) = u{r) is the interaction energy and U the potential of the box walls. 
The particles are polarizable and this leads to an attractive van der Waals interac- 
tion between the induced electric dipole moments. Such an interaction energy is 
proportional to r~ 6 , where r is the interatomic distance. At short distances, how- 
ever, the electrons repel because of the Pauli exclusion principle. Although the 
repulsive potential energy usually decays exponentially with r, we will write it in 
the form: 



u(r) 



oo 



-e(2r /r) e 



r < 2tq 
r > 2tq 



(2.5) 



(2r is the minimal distance to which the particles of "hard-core" radius r can 
approach.) 

The partition function of N interacting identical particles is 




Figure 2.1: Interatomic potential energy u(r) in real gases. 

Z(T,V,N) = j^J... J d* Pl ---d 3 p N J ••• J d 3 n---d 3 r N x 

e -^[E i (f?/ 2 ^+t / )+E i >,- . (2.6) 

After integrating over all momenta, the partition function becomes: 

Z(T, V,N) = ^ Sm Zl (T, V, AT), (2.7) 

where Zj carries all the interactions between the particles. We calculate the par- 
tition function in a mean-field approximation in which all the correlations are ne- 
glected. The potential acting on one particle is calculated by assuming a uniform 
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spatial distribution of other particles, 



Zi = 



fVi 

/ erne" 
Jo 



•/3$ 



N 



(2.8) 



V 1 is the volume, available to one particle, V x = V - (A^/2)(47r/3)(2r ) 3 = 
V — Nb. Here we subtracted from V the volume occupied by other particles. N/2 
is the number of pairs of particles. The mean-field potential $ is: 



$ = - / e (2r^rf(N/V)Anr 2 dr = 

J2r 



2Na 



(2.9) 



The resulting partition function is 

z(t,v,n) = ^( 2 ^ L ) m (V-Nbfc^ 



(a = (27r/3)e (2r ) 3 , b = (27r/3)(2r ) 3 ), the Helmholtz free energy is: 
F(T, V, N) = -k B T \nZ = -Nk B T In 



V -Nb f2nmk B T\ 3/2 



and the equation of state is 



p(T,V,N) = 



NknT N 2 a 



dV J T ,N V -Nb V 2 ' 
The last equation can be written in a more familiar form: 



(2.10) 



N 2 a 
V ' 
(2.11) 



(2.12) 



N 2 a 

(p+^)(V-Nb) = Nk B T. 



(2.13) 



This is the familiar van der Waals equation. We derived the van der Waals equa- 
tion of state in a mean-field approximation, we neglected all the correlations be- 
tween the molecules. Nevertheless, this equation qualitatively correctly describes 
the liquid-gas transition and it also predicts the existence of the critical point. 
However, as we neglected the correlations, which are essential in the vicinity of 
the critical point, it does not predict the correct critical behaviour. 
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The critical point is located at the inflexion point of the isotherm, (dp/dV) T = 
(d 2 p/dV 2 ) T = 0, and is located at: 

V c = SNb 
a 

Pc = m 2 

8a 

Below the critical point, the liquid and the gas phases are separated by a first- 
order phase transition. The equilibrium transition between the two phases takes 
place when their chemical potentials and therefore the Gibbs free energies are 
equal. 

Some comments are in place here. 

• Critical exponents. The van der Waals equation was derived in a mean-field 
approximation, therefore all the exponents are mean-field like. 

• Liquid — > gas transition. The liquid phase is stable until the point A is 
reached. In equilibrium, at this point, liquid starts to evaporate, the system 
enters the coexistence region (mixture of gas and liquid) until the point B is 
reached, which corresponds to the pure gas phase. 

If equilibrium is not reached, the liquid phase persists beyond the point A 
towards the point D. In this region, the liquid is metastable and can (at 
least in principle) persist until the point D is reached. Beyond this point, 
the compressibility is negative, the liquid state is unstable (not accessible). 

• Gas — > liquid transition. The situation is reversed. In equilibrium, the 
gas phase is stable until the point B is reached. This is the onset of the 
coexistence region which exists until the point A is reached. Metastable gas 
phase exists between the points B and E and, again, the region E to D is 
unstable. 



• Convexity. As we have discussed in Section 2.3.3 already, the Helmholtz 
free energy F of a stable system has to be a concave function of V. For 
metastable states it is concave locally whereas for unstable states, it is con- 
vex. 
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• Spinodals. The transition from a metastable (homogeneous) state toward 
the equilibrium (mixed) state is connected with the nucleation of droplets of 
the new phase. Droplets first nucleate and then grow until the equilibrium 
(coexistence) is reached between the new and the old phase. Nucleation is 
a dynamic process. 

The ultimate stability limit of the metastable state is called the spinodal line. 
Along this line the isothermal compressibility diverges, therefore the spin- 
odal line is in a sense a line of critical points. In the van der Waals (mean- 
field) systems the separation between the metastable and unstable states is 
sharp. Beyond MFA, for systems with short-range correlations, however, 
the transition between metastable and unstable states is often gradual. 

2.2 Ising Model 

The Ising model is the prototype model for all magnetic phase transitions and is 
probably the most studied model of statistical physics. In this model, the spin 

operator S is replaced by a number, which represents the z— component of the 
spin and is usually S — ±1 ("up" or "down"). The order parameter m — (S) is 
thus one -dimensional, it is a two-component scalar. 
The Ising Hamiltonian is 

n = -\ E JiASj - (2.i5) 

Jij is the interaction energy between the spins i and j. Usually (not always) the 
sum is only over nearest neighbouring (NN) pairs, then = J. H is the external 
magnetic field (in energy units). The case J > corresponds to ferromagnetism 
and J < to antiferromagnetism. As we shall show later, the Ising model has a 
phase transition at finite T if d > 1. In d = 1 the transition is at T = 0. Therefore 
we say that d — 1 is the lower critical dimension of the Ising model. 

In the low-temperature phase, the symmetry is broken, the spins are ordered, 
the order parameter is m ^ 0, whereas in the high-temperature phase the sym- 
metry is not broken, m vanishes, the spins are disordered. Therefore this model 
is a prototype of all order-disorder transitions. Other examples of order-disorder 
transitions, I would like to mention, are lattice gases and binary alloys. All these 
systems are equivalent to each other, they all have the same critical behaviour 
(which - of course - depends on d), for given d they all belong to the same univer- 
sality class. 
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2.2.1 Lattice Gas 

Imagine that a fluid (gas or liquid) system is divided into a regular lattice of cells 
of volume roughly equal to the particle (atom, molecule) volume. We say that 
the cell is occupied if a particle centre falls into this cell. Since the cell volume is 
comparable to the particle volume, there can be not more than one particle per cell. 
Notice that microscopic details (like positions of the particles) are irrelevant in the 
critical region, therefore the lattice gas model is suitable for studying the critical 
behaviour of real gases although the space is discretized. (It is also assumed that 
the kinetic energy of the gas particles, which is neglected in the lattice gas model, 
does not contribute a singular part to the partition function.) 

o o o m-m o o 

o o o o • o o 
o •-• o o o o 
o i-i o o • o 

o o o o o • o 

Figure 2.2: A configuration of a lattice gas in d = 2. Solid circles represent atoms 
and open circles empty sites. Lines show bonds between nearest-neighbouring 
atoms. 

A lattice gas is thus a collection of particles whose kinetic energy is neglected 
and that are arranged into discrete cells. The cells are either occupied by one atom 
[closed circles in Fig. 2.2] or empty [open circles]. The Hamiltonian is: 

Wig = J2 u ( r i,j)PiPj ~ vY^Pi, ( 2 - 16 ) 

where u is the interaction potential, 

r oo (nj = 0) 
u(rij) = I -u (r it j — NN distance) (2.17) 
^ (otherwise) 

and is the sum over pairs, pi is the occupation number, it is 1 if the site is 
occupied by an atom, and zero otherwise. The occupation is, like the Ising spin, 
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a double-valued function. Therefore we expect that the two models are related. 
Indeed, if we set 



the Hamiltonian 2.16 becomes: 



Here, z is the number of NN neighbours (4 for square lattice, 6 for simple cubic 
lattice), and Ho is a constant, independent of S. Thus, we see that Hi g is identical 
to the Ising Hamiltonian if we set u/A = J and (/i/2 + uz/A) = H. For u > 0, 
the atoms attract and they "condense" at low temperatures whereas for u < 
the atoms repel. At 1/2 coverage (one half of the lattice sites is occupied) and 
u < 0, the empty and occupied sites will alternate at low T, like spins in an 
antiferromagnet. 

The lattice gas model in d = 2 is used to investigate, e.g., hydrogen adsorbed 
on metal surfaces. The substrate provides discrete lattice sites on which hydrogen 
can be adsorbed. By varying the pressure, the coverage (amount of adsorbed 
hydrogen) is varied. This corresponds to varying the field in magnetic systems. 

There is an important difference between the lattice gas and the Ising models 
of magnetism. The number of atoms is constant (provided no atoms evaporate 
from or condense on the lattice), we say that the order parameter in the lattice gas 
models is conserved whereas it is not conserved in magnetic systems. This has 
important consequences for dynamics. 

2.2.2 Binary Alloys 

A binary alloy is composed of two metals, say A and B, arranged on a lattice. We 
introduce the occupations: 




(2.19) 




1 




if i is occupied by an A atom 
if % is occupied by a B atom 




1 




if % is occupied by a B atom 
if % is occupied by an A atom 



Pi + qi = 1. 



(2.20) 
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Depending on which atoms are nearest neighbours, there are three different NN 
interaction energies: uaa, ubb, and uab = uba- The Hamiltonian is: 

H ba = -U AA J2 PiPj - U BB J2 QiQj - U AB J2 (PiQj + 9*Pj)> ( 22V > 

ihj) (hj) 

where each sum runs over all NN pairs. If we write = (1 + Si)/ 2 and 
q i — {\ — Si)/2 then Si — 1 if % is occupied by A and Si = —1 for a B site. 
When there are 50% A atoms and 50% B atoms, J2i Si = 0, and the Hamiltonian 
becomes an Ising Hamiltonian in the absence of a magnetic field: 

7~iba = — J X! SiSj + ^o- (2.22) 

Ho is a constant, and J = (uaa + u B b — 2uab)/^- For binary alloys, the order 
parameter is also conserved. 

2.2.3 Ising model in d = 1 

It is not difficult to solve the Ising model in the absence of H in d = 1 exactly. 
The Hamiltonian of a linear chain of N spins with NN interactions is 

Hu = — J ^2 SiSi+i- (2.23) 

i 

We shall use periodic boundary conditions, that means that the spins will be 




Figure 2.3: Lattice sites of an Ising model in one dimension with periodic bound- 
ary conditions. 
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arranged on a ring, see Fig. 2.3. Then, the iV— th and the first spins are NN and 
the system is periodic. The partition function is: 

N 

Z = J2e- m{s}) = J2^ J ^ SiSi+1 = J2U^ JSiSi+1 - ( 2 - 24 ) 

{S} {S} {S}i=l 

(Don't forget: J2{s} * s me sum over a ^ s P m configurations - there are 2 N spin 
configurations - whereas J2i is the sum over all sitesl) 
The exponent in (2.24) can be written as: 

& PJSA +l = cosh p J SiS^+smhp J SiS^ = cosh (3 J+SiS i+1 sinh (3 J. (2.25) 

The last equality holds because SiSi + i can only be ±1 and cosh is an even function 
whereas sinh is an odd function of the argument. The partition function is: 

N 

Z = (cosh(3J) N ]T I] [1 + KSiS i+1 ] , (2.26) 

{S} i=l 

where K = tanh (3 J. We work out the product and sort terms in powers of K: 

Z = (cosh jf 2(1 + SiSzK^l + S 2 S 3 K) • • • (1 + SnS.K) 
{s} 

= (cosh £ [1 + ^(^1^2 + S 2 S 3 + --- + S N S!)+ 

{S} 

+K\S 1 S 2 S 2 S 3 + ...)+... + K N {S 1 S 2 S 2 ■ ■ ■ SnSnSJ] . (2.27) 

The terms, linear in K contain products of two different (neighbouring) spins, like 
SiSi+\. The sum over all spin configurations of this product vanish, Yl{ s} SiSi + \ = 
0, because there are two configurations with parallel spins (SiS i+1 = +1) and two 
with antiparallel spins (SiS i+ i = —1). Thus, the term linear in K vanishes after 
summation over all spin configurations. For the same reason also the sum over all 
spin configurations, which appear at the term proportional to K 2 , vanish. In order 
for a term to be different from zero, all the spins in the product must appear twice 
(then, J2{s t } &i = 2)- This condition is fulfilled only in the last term, which - af- 
ter summation over all spin configurations - gives 2 N K N . Therefore the partition 
function of the Ising model of a linear chain of iV spins is: 

Z = {2 cosh (3 J) N \l + K N ] . (2.28) 
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In the thermodynamic limit, N — > oo, K N — > 0, and the partition function be- 
comes: 



Z(H = 0) = (2cosh/3J) 



A? 



(2.29) 



The free energy, 



F/N 



Temperature 



Figure 2.4: Free energy of the Ising model in d = I. 

F = -Nk B T ln(2 cosh (3 J) , (2.30) 

is regular at all finite T, therefore the Ising model in d = 1 has no phase transition 
at any finite T. 

To investigate the possibility of a phase transition at T = 0, we will calculate 
the correlation function V and the correlation length £. We first replace J by Jj 
(at the end we put Jj = J), so that 

Z = 2 N JJcosh/3Jj. (2.31) 

j 

The nearest-neighbour spin-spin correlation function is: 

r(M + i) = (s t s t+1 ) = W^ s ^ E]JjSjS3+1 

Z {S} 

1 dZ 

The spin correlation function for arbitrary distance k, 

T(i,i + k) = (SiS i+k ) 
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— {SiSi+iSi+iSi^Si^Si+k-iSi+k-iSi+k) 
1 d d d d 

Zd(pJi) d(pj i+1 ) d{(3J i+2 ) ' ' ' 9(/3J i+fc _i; 



-Z, (2.33) 



r(fc) = r(z, z + fc) = [tanh(/3J)]' 



(2.34) 



is just the product of all the nearest-neighbour correlation functions between the 
sites % and i + k, and is (if all J's are equal) independent of the site %. The decay of 
the spin correlation function with distance is shown in Fig. 2.5. The correlations 
decay exponentially with distance and the correlation length is: 



1 



i = 

a lntanh/5J' 



(2.35) 



where a is the separation between the NN spins ("lattice constant"). In the limit 
T — > 0, the correlation length diverges exponentially as T — > 0: 



,2/3.7 



(2.36) 



a 2 

Equations 2.34 and 2.36 tell us that there is no long range order at any finite T and 



1 



2 4 6 8 10 
k 



Figure 2.5: Decay of the correlations in the d — 1 Ising model. 

that T = is the critical temperature of the d — 1 Ising model. Such a behaviour 
is typical for systems in the lower critical dimension: 

d > d c : T c ^ (2.37) 

d < d c : T c = 0, (2.38) 
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therefore we conclude that d — 1 is the lower critical dimension of the Ising 
model. 

The zero-field isothermal susceptibility of the d = 1 Ising model is, according 
to the fluctuation-dissipation theorem, 



XT(H = 0) = ^{S i S j ) 



(2.39) 



In the thermodynamic limit, xt becomes: 



Xt(H = 0) 



I 
N 

I 
N 



N/2 

i k=-N/2 
N/2 

i k=0 



(2.40) 



In the thermodynamic limit (N — > oo), Xt(H = 0) becomes 



XT (H = 0) = (3 



tanh (3 J 



^2J/k B T 



(2.41) 



The susceptibility has an essential singularity of the type exp(l/T) as T — > Tc = 
when J > 0. When J < 0, i.e., for antiferromagnetic NN interaction, the spin 
correlation function alternates in sign and the susceptibility is finite at all T. 
The internal energy is: 

E = -d In Z/df3 = -JN tanh f3J. (2.42) 

Notice the similarity between the internal energy and the short-range correlation 
function F(i, i + 1) ! This is a general property of the internal energy. 
The entropy is S = (E - F) /T, 

NJ 



S(H = 0) = Nk B In [2 cosh (3 J] - — tanh (3J. 



(2.43) 



In the limit T -> 0, 5 -> and in the limit T -> oo, 5 -> Nk B In 2, see Fig. 2.6. 
The heat capacity = T(dS/dT) H is: 

2 



83 



(2.44) 



^cosh /3J y 

The heat capacity is finite at all temperatures (no singularity at T = 0!), see Fig. 
2.7. 
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Figure 2.6: Temperature dependence of the entropy for the zero-field d — 1 Ising 
model. 




Figure 2.7: Temperature dependence of the heat capacity for the zero-field d — 1 
Ising model. 

2.2.4 Ising Model in d = 2 

The Ising model for H = in d = 2 dimensions has been solved exactly by 
Onsager. The calculation is lengthy, therefore we will first show that the Ising 
model in d — 2 has a phase transition at finite T and then quote the Onsager's 
result. 

Consider a two-dimensional square lattice with all the spins on the outer bound- 
ary "up". This simulates an arbitrarily weak external field and breaks the up-down 
symmetry of the model. To show that spontaneous ordering exists, it is enough 
to show that in the thermodynamic limit at a sufficiently small but finite tempera- 
ture there are more spins "up" than "down". Now we will estimate the fraction of 
"down" spins from the average size of a spin-down domain. 
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Figure 2.8: A possible domain configuration of the Ising model on a square lattice. 



With the above boundary conditions, all the domain walls must be closed 
loops, see Fig. 2.8. In a wall the neighbouring spins are antiparallel, this costs 
an energy 2 J per wall segment and the probability, that a domain with its total 
wall length (perimeter) b is thermally excited, is given by the Boltzmann factor 
exp[— 2(3 Jb]. For large b, there are many possible domain wall configurations 
with the same wall length b. For an estimate of the number m of possible wall 
configurations, assume that on each site the wall can either be straight or make 
a kink to either side. This makes a factor of 3 for each wall segment (except for 
the last one where it is bound to close the loop). Besides, a wall can start on each 
lattice site, this makes a factor of N. Therefore, m(b) is limited by: 



m(b) < N 



(2.45) 



and the average number of different thermally excited domains with perimeter b 
is: 



< 3 6 - 1 N e- 2/3Jb . 



(2.46) 



The number of "down" spins in a domain with perimeter b is iV_ < (b/A) 2 and 
the average concentration of all "down" spins is: 

N ~ ^ 16 

iV 6=4,6,8,- 1U 



1 d 2 

480a 2 



E 3 



b~—ab 



6=4,6,8,- 
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1 ^ 

1 d 2 ( q 2 \ 



48 9a x=2,3A 



48 <9a 2 VI -9/ 
q 2 (4-3q + q 2 ) 
12(1 -5)3 ' 



(2.47) 



where a = 2/3 J, x = 2b, q = 9exp(— 2a). For large, finite (3, q is small and 
the above ratio is < 1/2. This means that the order parameter (N+ — N-)/N is 
greater than zero at low T and the Ising model in d — 2 (and of course also in 
d > 2) has a phase transition at finite T. 

The exact critical temperature of the (anisotropic) Ising model on a rectangular 
lattice is given by the equation: 

sinh ( \ sinh ( \ = 1, 

\k B T c J \k B T c J 

where and J y are the spin interaction energies in x and y directions, respec- 
tively. 



2.3 Potts and Related Models 

In the previous Section we have discussed different classical two-state models 
(which are all mapped on the Ising model). A large class of classical physical 
systems cannot be described by the Ising model. As an example consider ad- 
sorption of noble gases on surfaces, like krypton on graphite. Graphite provides 
a two-dimensional hexagonal lattice of adsorption sites. Kr is large and once it 
occupies one hexagon, it is unfavourable for another Kr atom to occupy any of its 
NN hexagons, see Fig. 2.9. At full coverage, thus, Kr forms a regular lattice with 
1/3 of hexagons covered by Kr. Since there are three entirely equivalent positions 
(states!) for Kr (denoted by A, B, and C), the system is described by a three-state 
model. The adsorbed Kr forms A,B,C domains, separated by domain walls. It is 
clear from Fig. 2.9 that the A-B and the B-A walls are not equivalent, they have 
different energies in general. 

In the Potts model there are q states available to each site. The Potts Hamilto- 
nian (with only NN interactions) is: 

ft = -j£<W (2-48) 

(id) 
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Figure 2.9: Kr on graphite. There are three equivalent adsorption sites, but Kr 
atoms are too big to occupy all of them - only one third of the sites (i.e., either 
"A", "B", or "C" sites) can be occupied. We talk about "A," "B," or "C" domains. 
The domains are separated by domain walls. The A-B domain wall has a different 
structure (and energy) than the B-A wall. 



J is the interaction energy and 5 aijCrj is the 5-function which is = 1 if the NN sites 
% and j are in the same state and zero otherwise, 

_fla i = (T j with (ij) G NN 
\ otherwise K ^ 

Oi = 1, 2, • • • , q is a number which describes the state of the site i. Thus, in the 
Potts model, the interaction energy is equal to — J if the two NN sites are in the 
same state and vanishes otherwise. For q — 2, the Potts model is equivalent to the 
Ising model (S — ±1) because each site has two possible states in both cases. The 
two states can be visualized as two points on a line (which is a one-dimensional 
object). The mapping from the Potts to the Ising model is accomplished by writ- 
ing: 

Ss iSj = (2.50) 

For q = 3 the three states of the Potts model can be visualized as corners of 
a triangle (which is a two-dimensional object) whereas for q — 4, the four states 
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belong to corners of a tetrahedron, which is a 3-dimensional object. We conclude 
that the order parameter of the q— state Potts model is a q— component scalar. In 
d = 2, the Potts model has a continuous phase transition for q < 4 and a first-order 
transition for q > 4. 

Order parameter. Let x A , x B , • • ; x q be the fractions of cells (sites) that are 
in the states A, B, • • •, q, respectively. Of course, xa + xb + • • • + x q = 1. In 
the disordered phase, all the fractions are equal whereas in the ordered phase one 
state, say A, will have greater probability than the other, xa > (xb,xc, - ■ • ,x q ). 
We make a simplified approach, assuming that x B = x c = • • • = x q . The order 
parameter m measures the asymmetry in occupancies, 

Xb H hi, 

m = x A : 

q - 1 

1 + (q — \)m 1 — m 

XA = , X B = X C = ■ ■ ■ = X q = 

q q 

It is easy to verify that for rriA = 0: xa = xb = xc = • • • = x q and for m — 1: 
xa — 1 and x B = xc = ■ ■ ■ = x q = 0. thus, m A is the order parameter which is 
7^ in the ordered phase and vanishes in the disordered phase. 

I would like to mention another generalization of the 2-state model, that is the 
q— state (chiral) clock model. In this model, the states are imagined as hands of 
a clock (i.e., vectors on a circle) that can take only discrete orientations, 0j = 
(2%/q)ai (where <7j = 0, 1, 2, • • • , q — 1). The q— state chiral clock Hamiltonian is 

H = -JJ2 cos (©* - ©i + A • fij) (2.51) 

(iJ) 

A is the chirality which breaks the symmetry of the model, it produces an asym- 
metry in the interaction between the NN sites. The interaction energies between, 
say, <7j = 0, Gj = 1 and a,i = 1, aj = are different. The two-state clock model is 
identical to the Ising model and, if A = 0, the 3-state clock model is identical to 
the 3-state Potts model. For q > 3, the clock and the Potts model are not identical. 
The two-dimensional clock and Potts models are widely used in surface physics. 
The above mentioned Kr on graphite can be described with the 3-state chiral clock 
model with A^0. 

2.4 Two-dimensional xy Model 

A very interesting situation appears when the spin is a continuous vector in a plane 
(dimensionality of the order parameter is 2) and the spatial dimensionality is also 
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d = 2. As we have seen, the Ising model in d > 2 does have a phase transition 
at finite Tq- In the xy model, however, spin excitations (spin waves) are easier 
to excite thermally than for the Ising model, they are so strong that they even 
destroy long-range order at any finite T (Mermin- Wagner theorem). On the other 
hand, there has been evidence that some kind of a transition does take place at 
a finite Tc- What is the ordering and what is the nature of this transition? This 
question was answered in the early seventies by Berezinskii and by Kosterlitz and 
Thouless. 

The Hamiltonian of the xy model is: 

n = -Jj2^-S 3 = -Jj2 cos($, - $,), (2.52) 

where the spins are classical vectors of unit length aligned in the plane. $j is the 
angle between the spin direction and an arbitrarily chosen axis and is a continuous 
variable. For simplicity, let us consider the square lattice and let the sum be only 
over nearest neighbours. Only slowly varying configurations ($j — <C 1) give 
considerable contributions to the partition function, therefore we expand H: 

n-E = ^E(^-^) 2 = ^EIv$WI 2 

= J - J ' d 2 r\V${r)\ 2 . (2.53) 

Here a is the lattice constant and we transformed the last sum over all lattice sites 
into an integral over the area. $(r) is now a scalar field. 

Possible excitations above the (perfectly ordered) ground state are spin waves. 
They are responsible for destroying long-range order at any finite T. We will not 
consider spin waves here. Another possible excitation above the ground state is 
a vortex, shown in Fig. 2.10. Any line integral along a closed path around the 
centre of the vortex will give 



V$(7)-dr = 2nq .-.^ . (2.54) 




where q is the vorticity (q is integer; q — 1 for the vortex shown in Fig. 2.10) and 
4> is the polar angle to r. From the equations (2.53) and (2.54) we find the energy 
of this vortex: 

= f 2m *i(') , = * Jto (f)' > - ( " 5) 



E 
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Figure 2.10: Two examples of a vortex with vorticity q — 1. 

R is the size of the system. Notice that the energy of the vortex increases logarith- 
mically with the size of the system. In estimating the entropy associated with a 
vortex we realize that a vortex can be put on each lattice site therefore the number 
of configurations is about (R/a) 2 and the entropy of one vortex is 

S ^2k B \n(^j. (2.56) 

The free energy, F = E — TS, is positive and large at low temperatures. That 
means that it is not very easy to (thermally) excite single vortices. As temperature 
increases, the entropy term causes F to decrease. F vanishes (changes sign) at a 
critical temperature: 

k B T c = ^. (2.57) 

Above T c a large number of vortices is thermally excited. 

However, this does not mean that there are no vortices at low T. The easiest 
way to see this is to draw an analogy to a system of electrical charges in two 
dimensions. If k is a unit vector perpendicular to the plane of the system, then 
£ = V{q4>) x k is equal to the electric field produced by a charge of strength 
2-ntoq positioned at the centre of the vortex. The energy of an electric field is 

E el = € jJ d 2 r\£\ 2 = | / d 2 r|V(g0) x k\ 2 , (2.58) 
whereas the energy of a system of vortices is 

E vort = i [ d 2 r|V$ x k\ 2 = -E a . (2.59) 
Z J e 
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The energy of a system of vortices of strength ^ is equivalent to the energy of a 
system of electric charges! There is almost a complete analogy between a system 
of vortices and the two-dimensional Coulomb gas. We have seen that, like the 
electrostatic energy of a single charge in two dimensions, also the energy of a 
vortex diverges logarithmically with the size of the system. However, a pair of 
opposite charges, i.e., a dipole does have a finite electrostatic energy because at 
large distances the fields of two opposite charges almost cancel. In analogy to the 
energy of two charges ±q at fi and r 2 we write the total energy of a system of two 
opposite vortices: 

E 2 vort = 2irJq 2 In ~ ^ • (2-60) 

This self-energy of a vortex pair (of a dipole) is very small when the vortices are 
close together. Thus, it is difficult to excite a single vortex but it is very easy 
to excite a pair of opposite vortices. Whereas a single vortex influences the spin 
configuration on the whole lattice, the effect of a pair of opposite vortices cancels 
out on large distances from the pair so that the effect is limited to an area of the 
order \ fi — r 2 | 2 . 

Since the size of the system does not appear in (2.60), but it does apper in the 
expression for the entropy, we will have a finite density of bound vortex pairs. The 
equilibrium density of vortex pairs will depend on the pair-pair (dipole-dipole in 
the Coulomb gas case) interaction which we neglected here. 

Now we have the following picture. At low but > temperature (T < T c ) 
there will be a finite density of low-energy vortex pairs ("dipoles") of zero total 
vorticity. At high temperatures (T > T c ), on the other hand, there will be a large 
concentration of single vortices. The critical temperature, we estimated above, is 
thus associated with unbinding (dissociation) of vortex pairs. 

Very interesting is the behaviour of the spin correlation function (Si ■ Sj). 
We already mentioned that spin-wave excitations break the long-range order. If 
we neglect the spin wave excitations (which have nothing to do with the phase 
transition, they do not behave in a critical way at any finite temperature), then: 

i- /c c \ / > forT < T c n fin 

|^J^ ) = \o for T>T C . (m) 

This follows from the fact that bound vortices don't influence the spin orientations 
at large distances whereas they do if they are free. 

In a more realistic approach, one must take into account interaction between 
the vortex pairs, which was neglected in (2.60). Surrounding vortices (dipoles!) 
will screen the (electric) field and will thus reduce the self-energy of a vortex. 
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A sophisticated real-space renormalization group treatment of the two-dimensional 
xy model [Kosterlitz 1974] shows that the correlation length £ diverges above T c 
as: 

£ oc exp (b/t 1/2 ) 6wl.5 (2.62) 

where t — (T — T c )/T c . Below T c , however, the correlation length is always 
infinite and the correlations decay algebraically (like in a critical point) ! Therefore 
we can say that we have a line of critical points extending from zero to T c . The 
phase transition is characterized by a change from algebraic (power-law) decay of 
correlations below T c to an exponential decay above T c . 

Analogous to vortex unbinding transition in the xy model is the unbinding 
of opposite charges of a dipole in the two-dimensional Coulomb gas. At low 
temperatures, the Coulomb gas is composed of neutral dipoles, bound pairs, and 
is not conducting. At high temperature it is composed of free + and — charges 
which are conducting. 

Another analogy to the xy model are the dislocation lines in a (three-dimensional) 
crystal. At low temperatures, pairs of dislocation lines with opposite Burgers vec- 
tors can exist. At T c , the dislocations dissociate. This is the basis of the dislo- 
cation theory of melting. In this latter example, stiffness can play the role of the 
order parameter, it is ^ below the melting temperature and vanishes in the liquid 
phase. 

Let me stress again that the described behaviour of Coulomb gases or xy sys- 
tems is specific only for two-dimensions\ 

2.5 Spin Glasses 

The (vitreous) glasses are amorphous materials, their lattice is not crystalline but is 
strongly disordered. Topological defects destroy the long-range order in the crys- 
talline structure whereas the short-range order is still preserved to a high degree. 
Example: topological defects on a square lattice would be triangles or pentagons. 
Of course, the defects cause also a deformation of the squares. 

Spin glasses are the magnetic (spin) analogue of the (vitreous) glasses. In 
magnetic systems, the exchange interaction is often an oscillating function of dis- 
tance between the spins. If such a system is diluted with non-magnetic ions, the 
magnetic atoms are randomly positioned on the lattice and thus experience ran- 
dom exchange interactions with other, neighbouring atoms. Example: Aui_ x Fe x , 
where x is small. 
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Another system, which behaves in a similar way, is a random mixture of a 
ferromagnetic and an antiferromagnetic materials. The exchange interaction be- 
tween nearest neighbours % and j is random, it depends on the type of the magnetic 
atoms % and j. 

A Hamiltonian which reflects the essential features of a spin glass is: 

W = - E J(Ri,j)Si -Sj-^H- S,. (2.63) 

J(Ri,j) is a random variable and the first sum includes also more distant pairs than 
just nearest neighbours. Random exchange interactions cause that - even at low 
temperatures - the spins are not ordered, the system has no long-range order in 
the the ground state. The "usual" long-range order parameter 

i 

vanishes for a spin glass. Yet, in the experiment, a broad maximum in the specific 
heat and a rather sharp maximum in the zero-field susceptibility, connected with 
hysteresis and remanescence, are seen. So, there must be a kind of a transition 
as the temperature is varied, we shall call this a spin-glass transition and the low- 
temperature state the spin-glass phase. 

At high temperatures, the spins behave like in a normal paramagnet, they flip 
around dynamically so that the thermodynamic average (Si) vanishes. There is no 
long-range order in the system and also the local spontaneous magnetization, i.e., 
the (time) average of the spin at a site i vanishes, (Si) = for each site. At low 
temperatures, however, the spins freeze in a disordered configuration, the system 
is in a state where (Si) ^ although the average magnetization still vanishes, 
m = 0. The spin-glass transition is thus a freezing transition. 

Something happens to the dynamics of the spins, therefore we must consider 
time-dependent spins, S(t). To distinguish between the two phases, we introduce 
the (Edwards-Anderson) order parameter: 

q = lim lim 1 E(^o)^(t + *)>■ ( 2 -65) 

i 

q clearly vanishes in the paramagnetic phase and is > in the spin-glass phase, 
see Fig. 2.12. It measures the mean square local spontaneous magnetization. 
Notice that freezing is a gradual process, which starts at Tf where the first spins 
freeze and ends at T = where all the spins are frozen. 
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Figure 2.11: Spins in the low-temperature, spin-glass phase are frozen in a disor- 
dered configuration. 



Since, in the spin-glass phase the spins are (at least partially) frozen in a dis- 
ordered configuration, we can easily imagine that some spins are frustrated and 
that the spins do not freeze in a unique way. Freezing is history dependent, the 
configuration into which the spins freeze, depends on the conditions (e.g, on the 
magnitude of the magnetic filed) close to the freezing temperature. If we, upon 
cooling, apply an external field pointing "up", there will be more spins frozen in 
the "up" orientation and will stay in this orientation even if we switch the field 
off in the low-temperature phase. This means that the systems can fall into a 
metastable configuration from which it is not able to escape if the temperature is 
low. The free energy, thus, has many, many local minima in the configuration 
space, see Fig. 2.13. In the thermodynamic limit (N — > oo), some hills between 
the local minima will become infinitely high and the system will be trapped in 
such a minimum forever. In the time average, thus, the system will explore only a 
part of the total configuration space. This means that the time average of a phys- 
ical quantity is no longer equal to the ensemble average (where one explores the 
whole configuration space). A system in the spin-glass phase is thus not ergodic 
in the thermodynamic limit! Since it is ergodic in the paramagnetic phase, we say 
that at the freezing transition the ergodicity breaks down. 
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Figure 2.12: Temperature dependence of the order parameter q. Upon cooling, 
the spins start to freeze at Tf, below Tf, more and more spins are frozen until, 
eventually, at T = all the spins are frozen. 

2.5.1 Some Physical Properties of Spin Glasses 

The static susceptibility, Fig. 2.14, below the freezing temperature depends strongly 
on the way the experiment is performed. x(T) is largest and roughly temeprature 
independent in a "field cooling" (FC) cycle, i.e., if the field H is applied above Tf 
and then the sample is cooled at constant H below Tf. This susceptibility is re- 
versible, it is independent of the history. In contrast, the "zero-field cooled" (ZFC) 
susceptibility is obtained if the sample is cooled to T at H = and then the 
sample heated in the field. The difference between xfc and xzfc confirms the 
conjecture that the system can be frozen in different spin configurations (which, 
of course, also have different susceptibilities). 

The conjecture that freezing takes place at Tf is confirmed also by the ac 
susceptibility measurements, shown in Fig 2.15 for Cu 1 _ :E Mn :r (x = 0.9%) at dif- 
ferent frequencies. The lower the frequency, the more pronounced (less rounded) 
the cusp is. We see also that the position of the cusp shifts to lower temperatures 
as the frequency decreases. So, the freezing temperature should be defined as the 
maximum of x( v ) m me limit as v — > 0. 
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Figure 2.13: Free energy landscape in the configuration space of a spin glass at 
low temperatures. There is one global and many local minima. 




i 

15 

r(K) 



Figure 2.14: Temperature dependence of the static susceptibility of Cui^M^ for 
x — 1.08 and 2.02%. (b) and (d): cooled to low temperature at H = and then 
heated at H = 5.9 Oe; (a) and (c): cooled and/or heated in H = 5.9 Oe. [From 
Nagata et al., Phys. Rev. B 19, 1633 (1979).] 
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Temperature, T(K) 



Figure 2.15: Temperature dependence of the frequency-dependent susceptibility 
of Cui^Mn^ (x = 0.9%. [From C.A.M. Mulder et al., Phys. Rev. B 23, 1384 
(1982), ibid. 25,515 (1982).] 
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Chapter 3 
Quantum Models 



In classical statistical mechanics we did not know the exact microscopic state of 
the system (location in the phase space). We made use of the ergodic hypothesis 
and replaced the time average of a physical quantity with an ensemble average, 
i.e., an average over many equivalent systems. 

In quantum statistical mechanics the situation is similar. 



3.1 Quantum Statistics 

In this Chapter we will first introduce the density matrix and then make an intro- 
duction to the Bose condensation. 

In quantum mechanics, the energy levels are discrete and the state of the sys- 
tem is described by a wavefunction. Let us first assume that we know the dynami- 
cal state (wavefunction) \ip) of the whole system and the complete set of stationary 
states \n) which are the eigenstates of the full Hamiltonian of the compound sys- 
tem, 

U\n) = E n \n). (3.1) 

(We assume that the states \n) are normalized.) n denotes the quantum states (e.g., 
momenta) of all particles in the system. E n is the total energy of all particles of 
the system in the compound state n. Any dynamical state of the system can be 
expressed as a linear combination of these eigenstates, 

\m)=Y, c n(t)\n), (3.2) 

n 
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The coefficients c n (t) define a point in the (Hilbert) space of wavefunctions \n). 
The (time-dependent) wavefunction \tp(t)) which obeys the Schrodinger equation 

ih^- = H\1>) (3.3) 

is thus the quantum-mechanical analogue of the point in phase space of classical 
statistical mechanics. The expectation value of a physical quantity A (which does 
not necessarily commute with H) is: 

{A{t)) = (m\A\m) = £<c m (n|A|m> 

m,n 

(3.4) 

The last expression holds for any complete set of states \n), also if they are not 
the eigenstates of the full Hamiltonian. 



3.1.1 The Density Matrix 

As the number of particles in the system increases, the separation between the 
energy levels decreases rapidly, and, in the thermodynamic limit, it becomes ex- 
tremely high. On the other hand, the energy levels are never completely sharp. 
They are broadened for many reasons, including the uncertainty principle. In a 
macroscopic system (with extremely dense energy levels), the levels will thus al- 
ways overlap. Consequently, a macroscopic system will never be in a strictly sta- 
tionary quantum state, it will always be in a dynamical state, in a time-dependent 
mixture of stationary states. 

It is neither feasible nor possible to find a complete description of such a sys- 
tem. Like in classical statistical mechanics, we first assume that the system is 
ergodic, we replace the time average by an ensemble average over many quantum 
systems at the same instance. We talk about probabilities of finding the system 
in the states \n). The dynamical state is no longer a unique quantum state but 
is a statistical mixture of quantum states. We introduce the density matrix as the 
ensemble average of c* m c n : 

p = \m)p m;n {n\ 

p m ,n = (mlpln) = -Y,Cn*c l m . (3.5) 
" 1=1 

The sum over i runs over £1 configurations of the system, all being part of 
the same ensemble, and c l m = (m\ip l ). We do an ensemble average because our 
information on the system is not complete. 
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The density matrix is a statistical operator with the following important prop- 
erties: 

• The density matrix completely defines an ensemble of quantum systems 
and it carries all the information that is available for a quantum statistical 
ensemble. 

• Its diagonal elements p„ jn tell the probability that the system is in the state 
\n) whereas its non-diagonal elements p m n tell the probability of a transi- 
tion from the state \n) to the state \m). For stationary states, thus, p has to 
be diagonal, it commutes with the Hamiltonian. 

• The density matrix is normalized, 

Tr(p) = 1. (3.6) 

The trace runs over any complete set of states \n) and is independent of the 
choice of the basis set. 

• The (ensemble and quantum-mechanical) average value (A) of any physical 
quantity A is equal to: 

(A) = (ij\A\ij) = £ Pm,n(n\ A \™) = Tr(pA) = Tr(Ap). (3.7) 

n,m 

The trace in quantum statistical mechanics plays the role of integration over 
the whole available phase space in classical statistical mechanics. 

• The dynamics of the system is completely described by: 

ih^ = [H,p}. (3.8) 
This is the "Schrodinger equation" of the density matrix. 

3.1.2 Ensembles in Quantum Statistical Mechanics 

Now, we shall specify the density matrix and its relation to thermodynamic quan- 
tities. 
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Quantum Microcanonical ensemble. 

The energy of systems belonging to a microcanonical ensemble is fixed, 
therefore it is convenient to use the energy eigenstates as a basis. Then the 
density matrix is diagonal. According to the postulate of equal weights, 
all the states \n) with the energy E n between E and E + AE are equally 
probable and let there be AVI such states. For a quantum system, Ail is an 
integer number. The probability for the system being in one of these states 
is p = 1/Afi and the entropy is equal to: 



S = -k B (\np). 



(3.9) 



In general, we must use the density matrix instead of simple probabilities p 
and the entropy of a quantum system is: 



S = —k B (\np) = —ks Tr(plnp) 



(3.10) 



Again, the average ( ) means ensemble and quantum average - both aver- 
ages cannot be separated. 

Canonical ensemble. In the energy representation p is diagonal and we can 
use the same arguments as in the case of classical canonical distributions. 
The matrix elements are equal to: 

1 



Pr, 



-PE„ 



where the normalization constant is (again) the partition function 



-PE„ 



(3.11) 



(3.12) 



In an arbitrary basis, the density matrix can be written more generally as 




with 



Tre 



-PH 



(3.13) 



(p and H are operators!). Thus, for a canonical ensemble, the average (ther- 
mal and quantum- mechanical) of an operator A is: 



1 



(A) = Tr (pA) = - Tr (e"^ A) 



(3.14) 
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In particular, the internal energy is 



E=(H) = l-Tx (He-*™) = - A i nT r e"^ 
/' ^ up 



The free energy is 



E=-^- j \nZ. (3.15) 
op 



-k B T\n Z. (3.16) 



All the thermodynamic relations are the same as before (Tables I - III are 
still valid), the only difference is that now we have to calculate the partition 
function as the trace of the density matrix. 

The entropy in a canonical ensemble, S, is the ensemble average over all 
(energy) states of the microcanonical entropy S: 

S = -(k B \n p) = -fc B Tr(plnp). (3.17) 

• Grand canonical ensemble. In the grand canonical ensemble the operator p 
operates on a generalized Hilbert space which is the direct sum of all Hilbert 
spaces with fixed number of particles. The density matrix is: 



where the grand partition function is 

~(T,V,p) =Tre-^ n ^ N \ 

Notice that now iV is an operator in the generalized Hilbert space. As in the 
classical case, S is equal to the canonical partition function Z, weighted by 
the fugacity z = exp((3p) and summed over N. The ensemble average of 
an operator A in the grand canonical ensemble is 
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3.2 Bose Systems 



3.2.1 Ideal Bose Gas 

Now we shall consider the properties of an ideal gas of particles with integer spin 
and mass m > 0. As usually for ideal gases, we neglect the interaction between 
the particles. The Hamiltonian of the ideal Bose gas is: 



H 



^ 2m v 

v 



np is the number of particles in the state with momentum p. Because we are deal- 
ing with Bosons, can also be > 1. We want to study condensation where two 
phases come into contact, the number of particles in one phase is not fixed, there- 
fore we must work with the grand-canonical distribution. The grand-canonical 
partition function of ideal Bose gas is: 



(Here we will disregard the (2S + 1) factor which comes from the spin degener- 
acy.) In the grand partition function, each sum over is over all non-negative 
integers and, since the particles do not interact, we can split the partition function 
into a product: 



s= E E 

no=0 ni=0 



no 



-/3(pl/2m-a) 



(3(p 2 2 /2m-fM) 



"2 



■}= (3.18) 




-/3(p 2 /2m- fj,) 



n- 



ze -/3p2/2m ' 



Here we have introduced the fugacity z = e^. The equation of state of an ideal 
Bose gas is (see Table I): 

PV = k B T\nE(T,V,fi) 



PV 



-J> [l-C 



-/3(p 2 /2m-M) 



(3.19) 
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The thermal average of N is determined by the chemical potential // (see Tables): 
We can write (JV) = J2pnp where 

= e /3(p2/2m- M ) _ X ( 3 - 21) 

is the Bose-Einstein distribution function, which tells the occupation probability 
of a given (non-degenerate) state p. For (N) to be finite and positive, the second 
term in the denominator has to be < 1 for any p, in particular for p = 0. This 
means that /j, must be < and < z < 1. 

For large V we replace the sum by an integral: 



and we get: 



47rl/ /*°° 2 , 1 

= ^3/2^) (™^§f) = J #3/2(2), (3-23) 



where 



1 /-oo r ^-l J 

= ^_ / dx— — = Y - (3.24) 

y K ' r(v) Jo ^z- 1 - 1 f-r Z" 



i» y 

and A is the thermal wavelength, 

/ 27Tft 2 \ V2 

A 



rnksT 



The functions # 3 / 2 (2) and #5/2(2) are shown in Fig 3.1. For later use: #3/2(1) 
2.612 and #5/2(1) 1.342. The pressure is: 

' /,; ^ (A; B T) 5 / 3 # 5 / 2 (^) = -^#5/2(2). (3.25) 



V2tt^V v " 7 y&//v 7 A3 
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Figure 3.1: The functions g 3 / 2 (z) and g 5 / 2 (z). 



Now let us start increasing the pressure at fixed (high enough) T. This increases 
the chemical potential, the fugacity, and the density of particles, (JV) /V. How- 
ever, we know that z must be < 1. What happens when z — > 1? Can we further 
increase the pressure or the density of particles? 

To answer this question we must go back to Eq. (3.20). The term with p = 
which was not included in (3.23) (The density of states n(E) has a \[E depen- 
dence on E and vanishes as E — > 0.) becomes singular (divergent) when p — > 
and it has to be treated separately. This divergence has very important conse- 
quences for a Bose gas, as we shall show now. Instead of Eq. (3.20) we must 
write: 



<ao = E 



1 



m e 



f3(p 2 /2m-n) 



V 



(3.26) 



and the equation of state is 

PV AttV 



k B T 



h 3 



p 2 dp In 



1 - ZQ -^/2m 



-ln(l -z) 




(3.27) 
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The last terms in Eqs. (3.26) and (3.27) come from the termp = and correspond 
to the particles in the lowest energy, in the ground state. In the thermodynamic 
limit the last term in (3.27) vanishes whereas it remains finite in (3.26). We will 
denote 



(N) = (N') + (N ), (3.28) 

where (N 1 ) is the average number of bosons in the continuum (excited, p > 0) 
states, 

(N') = y 3 g 3/2 (z) (3.29) 

and {N ) the average number of bosons in the ground state, 

(No) = (3.30) 

The average density of bosons in the continuum states reaches its maximum value 
when gs/2(z) is maximal, that is for z — > 1 Qx — > 0), 

^%^= 93 , 2 (l)(ff) 3/! «T^. (3.31) 

The fact that (N') is limited and cannot exceed (N' max ) has its origin in the 
quantum-mechanical and statistical nature of the system. The number of avail- 
able states in a box with volume V is finite and is equal to j^. Each state, on the 
average, accommodates n$ particles, where is given by the Bose-Einstein dis- 
tribution function, Eq. (3.21). If there are more than (N 1 ) particles in the system, 
they are pushed into the ground state. 



3.2.2 Bose Einstein Condensation 

We will now investigate the behaviour of ideal Bose gases, described by Eqs. 
(3.26, 3.27). Let us consider a system of N bosons in a volume V. 

At high enough temperature, (N' max ) is so large that iV < (N' max ) and the 
chemical potential is determined from 

N =^= V ( ! t§-f' 1 9sdz). (3.32) 
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This means that z < 1 and // < 0. The density of particles in the ground state is 

(N )/V = */(l -z)/V -> as \/^oo. (3.33) 

Almost all the particles are in the |pj > states and the ground state is macro- 
scopically empty. The equation of state simplifies to 

P _*/>(*) (334) 



k B T A 3 

At very high T, g 3 / 2 (z) and #5/2(2) are 1 and we approximate [see the series 
for n ,Eq. (3.24)]: 

#5/2(2) « 5-3/2 « 2 (3.35) 
and the equation of state becomes: 

P .Jmk B T\ m (mk B T\ m z 

z is determined by the density of particles: 

N (mk B T\ 3/2 (mk B T\ 3/2 z 

The fugacity 2 = (N)\ 3 /V vanishes as T -3 / 2 at high T and for fixed density. 
This justifies the approximation (3.35). Eliminating z from the last two equations 
yields the equation of state of the ideal Bose gas at high T: 



PV =< N > k R T 



(3.38) 



which is identical to the equation of state of a classical ideal gas (as it should be 
at high T!). 

Now we lower the temperature by keeping N and V fixed. Eq. (3.32) tells 
us that 53/2(2) must increase. This means that z and fi also increase. Eventually, 
a temperature is reached, where z reaches its maximum value (z — > 1, N = 
(N' max )). This is the transition temperature T c and is determined by the condition 

N f mk B T c \ 3/2 

v = */ 2(1 H~2^H (339) 
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from which we find: 



2nh 2 



k B m[g 3/2 (l)}^ \ V 



A\ 2 / 3 



(3.40) 



If we insert this expression into the equation of state, Eq. (3.34), we obtain the 
pressure at which the transition takes place: 



P c = k B T c /2 g 5/2 (l) 



mk 



B 



3/2 



oc T, 



5/2 



C 



(3.41) 



At Tc, the ensemble average of the density of particles in the continuum states 
reaches its maximum value whereas the ground state is still empty, (N )/V = 0. 

Below T c , (N' max ) further decreases and becomes < A (we keep A fixed). 
The chemical potential further approaches zero (/x cannot be exactly = 0, because 
that would mean that A = oo) and z — > 1 in such a way that the total number of 
particles is kept constant, z is very close to I, we can put g n (z) as g n (l) and the 
number of particles in the continuum states is 



(N 1 ) = V 



' inksT 



3/2 



^3/2(1)- 



(3.42) 



When (N') < N, the ground state starts to fill, (jV ) increases. Below T c , thus, 
only a part of the particles can be accommodated in the continuum states. The rest 
must go into the ground state! The number of bosons in the ground state is: 



(N ) — N — (N 1 ) = N - V 



rnksT 



3/2 



53/2(1) 



from which we get: 



N 



Tc) 



3/2 



(3.43) 



(A^o) /A is finite, the state with p = is occupied with a macroscopic number of 
particles, the particles condense in the momentum space into the zero-momentum 
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Figure 3.2: Isotherms of the ideal Bose gas. 




0.5 1 1.5 

T/Tc 



Figure 3.3: Temperature dependence of the relative number of bosons in the 
ground ((N ) /N) and excited ((TV') /N) states. 

state. This is the Bose-Einstein condensation. The fugacity is, from Eq. (3.30): 
z = (N ) /((N ) + 1). In the thermodynamic limit, when V and (N ) — > oo, 
z = 1 and /i = below T c . The equation of state is: 

h Jmk B T\ 3/ \ g 5/2 (l) (N f ) 

Only the particles in the continuum states ( in the gas phase ) contribute to pres- 
sure. The particles in the ground state (condensate) are at rest, they cannot exert 
any pressure. If the density is increased, the extra particles fall into the ground 
state and the pressure does not increase. At fixed T <T C the density of particles 
in the gas phase is constant and P is independent of V, see Fig. 3.2. The temper- 
ature dependence of (No) is shown in Fig. 3.3. In the limit T — > 0, (N' max ) — > 
and all the particles are in the ground state, N = (No). 
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Gas phase 
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Figure 3.4: Phase diagram of an ideal Bose gas. 



The corresponding phase diagram is shown in Fig. 3.4. The normal, gas phase 
exists for T > T c , i.e., to the right of the transition line. On the line, given by 
(3.41), condensation takes place and since the particles in the ground state don't 
contribute to pressure, the condensate lies on the transition line itself. Notice that 
the "critical points" in Bose systems form a line in the P vs. T or in the P vs. V 
planes and are not in a single point like in van der Waals gases. 

We now invert Eq. (3.44), 



As T — > 0, also (N') — > and therefore V — > at any finite pressure. The 
condensed phase does not occupy any volume. This means that the ideal Bose gas 
can be compressed to zero volume without any increase in pressure. 

Some comments: The condensate (ground state) contributes neither to E, nor 
to Cy, P, or VI At low T, C v vanishes like T 3 / 2 in contrast to the photon gas, 
where Cy oc T 3 . They differ because the density of states is different in both 
cases. There are more excited states available for particles than for photons and 
the specific heat is greater. 

Of course, this, and infinite compressibility are the artifacts of the ideal gas 
model. In reality, because of atomic repulsion, the volume of the condensed phase 
does not vanish and the compressibility does not diverge. The equation of state of 
a non-ideal Bose gas is shown in Fig. 3.5. 



V = 



05/2(1) (N') 
93/2(1) P 



k B T. 



(3.45) 
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Figure 3.5: An isotherm of a non-ideal Bose gas. 




point 
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Figure 3.6: Phase diagram of 4 He. He I is the normal, and He II the superfluid 
phase. 



C0 A - 




Figure 3.7: Dispersion curve of elementary excitations in liquid 4 He. The excita- 
tions around the minimum at p are called rotons. 
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3.3 Liquid 4 He 



A candidate for the Bose-Einstein condensation at low T is liquid 4 He. The phase 
diagram is shown in Fig. 3.6. He I is the normal liquid and He II exhibits super- 
fluid behaviour. From the fact that the line separating the solid and the fluid phases 
is horizontal at T — > 0, we conclude with the help of the Clausius-Clapeyron 
equation, 

dp „ •S'fluid — SUiki ,~ At ~. 

^r =0= t7 t7 — (346) 

al ^fluid — ^solid 

that liquid He II at T — > has no more entropy than solid helium, which is highly 
ordered, of course. We conclude that the particles in the ground state (T — > 0) 
have already the highest possible degree of order although they are in the fluid 
and not in the solid phase. We interpret the He II phase as a mixture of two fluids, 
a normal fluid with final viscosity and a super-fluid (formed of particles in the 
ground state) with zero viscosity. Characteristic for the superfluid state is also a 
very large heat conductivity (as a consequence, He II does not boil, the atoms only 
evaporate from the surface). What is the ground state of He II and what are the 
excitations? 

We shall first discuss possible excitations above the ground state of liquid 4 He 
and only then the properties of the ground state. In the liquid, the atoms are 
close together, they are strongly interacting (exchanging momenta). A state of 
an atom with fixed momentum is not an eigenstate of the Hamiltonian, the atoms 
are subject to a collective motion. The elementary excitations in the system are 
collective excitations, quasiparticles, phonons (longitudinal because it is a liquid) 
and rotons (see Fig. 3.7). At low momenta, the excitations are phonons with a 
linear dispersion relation, 

e p = cp. (3.47) 

At high momenta (short wavelengths), the interaction between the atoms causes 
a minimum in the dispersion curve. Elementary excitations around this minimum 
are (for historical reasons) called rotons. Their energy is written in the form: 



(P ~ Po) 2 
2m* 



e p = A+ yF „_ F y , (3.48) 



where A is the energy gap, A/k B = 8.6K, and m* = 1.0 x lCT 27 kg is the 
effective mass of rotons. po = 2.0 x 10 _24 kg m/s. At low T, the density of these 
quasiparticles is low and they don't interact. 
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The grand canonical partition function for phonons (p = 0) is, in analogy 
with (3.18): 

Ins,* = J o p 2 dpln(l-e-^) 

AttV r°° x 3 dx 



3(h(3c) 3 Jo e x - 1 

- |S> T > 3 (3 - 49 > 

The grand canonical partition function for rotons (also p = 0) is: 

47rV r 00 



lnS rot - - — 



/ 00 j9 2 ^ln(l-e-' 3A - /3 ( p - p «) 2 / 2m *) 



At low T (large /?), e ^ A ^ p po) 2 / 2 " 1 * is small and we approximate ln(l — x) 
—x. (This is equivalent to using the Maxwell-Boltzmann distribution.) 



l nErot = ^ e ^ A r P 2 dpe-^ p - p0)2/2m * 
h 6 Jo 

/oo 



— e^ A /°° v 2 dve- p{ *- po)2/2 "> 



= "W~ ( P » + Tj e (3 ' 50) 

From these equations, one can calculate the specific heat (and other thermody- 
namic quantities) at low temperatures. 



3.3.1 Superfluidity 

Consider an isolated system (in which energy and momentum are conserved) com- 
posed of a long cylindrical capillary filled with liquid 4 He at T = 0. Initially, let 
the capillary tube rotate along the symmetry axis with the tangential velocity v 
and let the fluid be at rest. We would expect that - because of friction - energy 
would be transferred from the tube to the fluid and that the fluid will gradually be 
excited out of the ground state. Let us check under which conditions the energy 
and momentum can be transferred to the fluid. Immediately after the tube started 
to rotate, the tube had the kinetic energy (M is the mass of the tube) and 

the fluid was at rest and had the ground-state energy E . At a later time, we expect 
the tube velocity to decrease by Av (because we consider an isolated system) and 
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that in the fluid an (at least one) elementary excitation with energy e p and momen- 
tum p is created. The conservation laws demand that the energy at a later time 
is: 

-M(v - Av) 2 + e p + E = -Mvl + E o 

(3.51) 

and that the momentum is: 

M(v - Av) +p t = Mv =>• MAv = p t (3.52) 

where p t is the excitation momentum in the tangential direction. For Av ^ 0, the 
first equation gives the inequality 

e p = Mv Av - \m{Av) 2 < Mv Av. (3.53) 

After eliminating MAv, the conservation laws tell us that 

v-p>e p . (3.54) 

Inspection of Fig. 3.7 shows that this is only possible if v is larger than a threshold 
velocity v c . At small relative velocity between the tube and the fluid, the rotons 
are not excited, the momentum is not transferred to the fluid, the fluid will stay at 
rest. This means that there is no friction, no viscosity - this is the superfluid phase. 
(As we will see later, roton excitations are not responsible for the excitation out 
of the ground state. Vortices have higher energy, but lower e/p ratio.) 



3.3.2 Vortex Excitations 

There are other collective excitations in liquid 4 He that have higher energies than 
phonons or rotons, but which nonetheless have much lower energy to momentum 
ratio than rotons. The creation of these excitations allows energy and momentum 
to be transferred to the fluid at a much lower critical velocity than in the case of 
rotons. 

Let us now discuss such excitations of the condensate. The wave function of 
an excited collective state is 

* = ^ (r) , (3.55) 

the (mass) current density is 

-> ih 



hi/j 2 V(j)(r). (3.56) 
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On the other hand, 



A 7 



(3.57) 



3 = 



V 



mv 



where N/V = p = ip 2 . From the two equations we get for the particle velocity 

i?=— V0(r). (3.58) 
m 

The wavefunction must be single-valued and the phase is determined modulo 
2%, therefore the closed-path integral of the momentum is 



where h is the Planck's constant and n an arbitrary integer. 

To understand the nature of these excited states it is best to consider liquid 4 He 
within a small rotating bucket. If the bucket is slowly brought into rotation, the 
liquid will behave in the following way. At small angular velocities of the bucket, 
the liquid will remain at rest. If the angular velocity is increased beyond some 
threshold angular velocity, the energy is transferred to the fluid, and the fluid will 
make an abrupt transition to a state in which there is a vortex line running from 
top to bottom, parallel to the axis of rotation. The vortex line looks like a small 
whirlpool, the fluid flows around the line with the velocity that decreases with 
distance from the line. The vorticity, defined as the line integral of the momentum 
along any path encircling the vortex line, Eq. (3.59), is quantized. As the angular 
velocity is increased, more and more quantized vortex lines are created. Once a 
vortex line is created, the macroscopically occupied single-particle ground state 
is no longer a zero-momentum eigenstate. It is a quantum state with a flow pat- 
tern that includes the vortex lines. When there are many vortex lines, they will 
distribute themselves in a pattern so that the average fluid velocity (averaged over 
many vortices) is the same as in the case of rigid-body rotation. 

To summarize, the quantization of vortices is related to the phase factor of the 
condensate wave function, defined in Eq. (3.60). In the lowest-energy state, ip and 
do not depend on f, they are constant. The state in which varies with r is an 
excited state which describes the fluid with one or several vortices. 

3.3.3 Order Parameter 

A superfluid can be considered as a mixture of two coexisting fluids, of a normal 
component and of a superfluid component. At a first glance, one would choose 




(3.59) 
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the density of particles in the ground (superfluid) state as the order parameter, 
because it is 7^ in the low-temperature phase and vanishes at Tq. However, in an 
interacting system, the operator N doesn't commute with H, the single-particle 
momenta are not good quantum numbers. On the other hand we already know 

— * — * 

that the ground state is a state with zero momentum p = = — ihV V = 0. 
This means that the wavefunction is a constant in the ground state. Therefore we 
choose for the order parameter the statistical average of the ground-state wave 
function which commutes with H and which can be written in the form: 

(tf ) = - Tr ^jQ-m-^ = ^ (3 _ 60) 

(ip is the amplitude and the phase of the wavefunction). The density of the 
superfluid component is (N )/V = i\) 2 = constant. In the superfluid phase, a 
macroscopically occupied single -particle quantum-mechanical ground state with 
zero momentum extends through the whole volume available to the fluid. Each 
atom is not at any particular place, it is simultaneously everywhere in the macro- 
scopic volume! This is the meaning of the order parameter defined in (3.60). 
Symmetry breaking. There is a strong parallel between (\&) and the sponta- 

— * 

neous magnetization (M) in ferromagnets, 

(M) = Tr (Me _/W ) , (3.61) 

Zj ^ ' 

— * 

where M is the (operator for the) total magnetic moment and 

H = —J Si ■ Sj. 

In the absence of an external magnetic field, thus, H is invariant under rotations, 
M and — M occur with equal probability and therefore (M) is always zero. But 
we know that the ferromagnets do have spontaneous magnetization at low T. The 
resolution lies in "spontaneous symmetry breaking" - the ground state of a Hamil- 
tonian does not possess the full symmetry of the Hamiltonian. The ground state 
wavefunction is not rotationally invariant, because (M) points along a definite di- 
rection in space whereas H is rotationally invariant. In the case of ferromagnets, 
thus, the transition corresponds to spontaneous breaking of the rotational symme- 
try. One must therefore redefine (M) in such a way that +M and — M are not 
both included in the partition function with equal weights. This is most easily 
done by imposing an infinitesimal external field H: 

(M)= lim lim -TrMe - ^-^. (3.62) 

H^O V-+00 Z 
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Returning to Bose systems, the Hamiltonian has a global gauge invariance, 
the wave functions are invariant with respect to the phase shift, 



(3.63) 



where is an arbitrary phase factor (real number). Bose condensation corresponds 
to a spontaneous breaking of global gauge invariance. When calculating the order 
parameter, we imagine, like in magnetic systems that the system is subjected to an 
(infinitesimal) external field rj which couples linearly to the order parameter. We 
calculate the statistical average of ^(f) in the thermodynamic limit: 



(*(r)> 



lim lim 



1 



Tr*(r)e' 



-0£ 



(3.64) 



where 



£ = H- fiN - I d 3 r[ty(r)ri(f) + &(r)rjt(r)]. (3.65) 

The only essential difference with the ferromagnetic case is that the external field 
here is fictitious, it cannot be realized experimentally. 
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4.1 Mean-Field Theories 

Very few models of statistical mechanics have been solved exactly; in most of 
the cases one has to rely on approximative methods. Among them, the mean-field 
approximation (MFA) is one of the most widely used. The advantage of the mean- 
field theory is its simplicity and that it correctly predicts the qualitative features of 
a system in most cases. 

The essence of the mean-field theory is the assumption of statistical indepen- 
dence of the local ordering (spins in the case of magnetic systems). The interac- 
tion terms in the Hamiltonian are replaced by an effective, "mean field" term. In 
this way, all the information on correlations in the fluctuations is lost. Therefore, 
the mean-field theory is usually inadequate in the critical regime. It usually gives 
wrong critical exponents, in particular at low spatial dimensions when the number 
of NN is small. 

The MFA becomes exact in the limit as the number of interacting neighbours 
z — > oo. This is the case when the range of interaction z — > oo or if the spatial 
dimension is high. In both cases, a site "feels" contributions from many neigh- 
bours and therefore the fluctuations average out or become even irrelevant (for 
Ising systems, they become irrelevant in d > 4). At low d, the MFA must be used 
with caution. Not only that it predicts wrong critical exponents, its predictions are 
even qualitatively wrong sometimes. For example, it predicts long-range order 
and finite critical temperature for the d = 1 Ising model. In the following, we 
shall introduce and formulate the MFA for magnetic systems. In the Appendix A 
we summarize basic thermodynamic relations for magnetic systems. The results, 



73 



I. Vilfan 



Statistical Mechanics 



Chapter 4 



however, are valid also for other systems. 

Two formulations of the mean-field approximation. We will now formulate 
the mean-field approximation in two different ways. 

• A: The "Weiss molecular field theory" according to its "inventor" Pierre 
Weiss. This method is straightforward and easy to understand. It gives us 
an expression for the order parameter but does not tell us anything about the 
free energy or partition function. 

• B: The second method will be the Bragg-Williams approximation which is 
more sophisticated and is based on the free energy minimization. 

We will develop the mean-field approximation on the Ising model in an exter- 
nal magnetic field. 

4.1.1 Weiss Molecular Field of an Ising System 

Weiss formulated a theory of ferromagnetism in which he assumed that the effect 
of the neighbouring spins on a given spin can be described by a fictitious "molec- 
ular field" which is proportional to the average magnetization of the system. 
We start with the Ising Hamiltonian in an external magnetic field: 

(ij) i 

The local field acting on the spin at the site % is: 

Hi = Jj2Sj + H - (4-2) 

j 

It depends on the configuration of all the neighbouring spins around the site i. 
If the lattice is such that each spin has many neighbours, then it is not a bad 
approximation to replace the actual value of the neighbouring spins by the mean 
value of all spins, 

Y^S^zm, (4.3) 

3 

where z is the number of nearest neighbours and m is the mean (average) magne- 
tization per site, m = ~^J2j Sj- With (4.3), the local field becomes 

Hi = zJm + H, (4.4) 
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Figure 4.1: Temperature dependence of the order parameter in the mean-field 
approximation. 

The actual local field acting on the site % was replaced by its mean value and is 
independent of the site, Hi — > H M f- In the mean-field approximation the original 
Hamiltonian is replaced by a mean-field Hamiltonian 



which is equivalent to the Hamiltonian of a non-interacting spin system in an 
external magnetic field and can be solved exactly: 



This is a self-consistent, transcendental equation for m. It has a non-trivial solu- 
tion (m 7^ 0) for (3z J > 1, i.e., at low temperatures. The temperature dependence 
of the order parameter m for H = is shown in Fig. 4. 1. It vanishes at the critical 
temperature Tc = zJ. 

4.1.2 Bragg- Williams Approximation 

The Gibbs free energy is a function of T and H, 



7~Lmf — — 




(4.5) 



m = tanh(/5 H M f) = tanh[/3(z Jm + H)]. 



(4.6) 



dG = -SdT - MdH. 



(4.7) 
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This equation tells us that for fixed T and H, 5G = 0. In equilibrium, thus, G 
is an extremum of a trial free energy Q with respect to any variational parameter. 
We will now search the trial free energy Q. 

We start with the density matrix and calculate the Gibbs free energy of an 
Ising system. In the Bragg- Williams (mean-field) approach we approximate the 
total density matrix p by a direct product of the density matrices for individual 
spins, 

P^UpitxUe-^. (4.8) 

i i 

Physically this means that the spin fluctuations are considered to be uncorrelated, 
the orientation of a spin is statistically independent of the orientation of any other 
spin. There are no correlations in the spin orientation except those coming from 
the long-range order. Under this assumption, the Ising Hamiltonian is not only 
diagonal in the spin space (the state Si = +1 is uncoupled from the state Si = 
— 1; physically this means that the Ising system has no dynamics), but it is also 
independent of all other spins in the system. Thus, the trial density matrices p, 
are diagonal and of the general form (remember, its trace must be = 1): 



Ul + m') 



2 



\{l -m'i) 



(4.9) 



They depend on the variational parameter which is site-independent, m! i = ml, 
because of translational invariance. With this density matrix the average spin 
orientation is: 

{S i )=Tr(S iPi )=m. (4.10) 

Because of statistical independence of the spins in MFA, the short-range spin 
correlations are: 

(S i S j ) = (S i )(S j )=m 2 . (4.11) 
We can write the (trial) enthalpy E = (H) [where H is given by the Eq. (4.1)] as: 

E = -^zJNm' 2 - NHm' (4.12) 

and the trial entropy as: 

, m t i \ at, fl + m' 1 + m' l—m' 1 - m'\ 
S = -k B Tr(p In p) = -Nk B — - — In — 1 — In — - — 
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Figure 4.2: Trial Gibbs free energy g(T,H;m') vs. w! at H = 0.1 z J. Solid 
circles denote stable states, open circle a metastable state, and the solid triangle an 
unstable state. The curves are labelled by T/T c . m and G(T, H) are determined 
by the minimum of Q. 

(4.13) 

With these expressions, the trial Gibbs free energy becomes: 

g(H,T;m') = E-TS = E + k B T Tr(plnp) = -\zJNm' 2 - NHm'+ 

V m' l + m' 1 — m' 1 — m' \ A A 

In 1 In (4.14) 

2 2 2 2 ; v ' 

It depends on H, T, and on the variational parameter ml. Q for H = 0.1 (in 
units of z J) is shown in Fig. 4.2. To get the true equilibrium Gibbs free energy, 
G(H, T; m) has to be minimized with respect to the variational parameter ml, 

G(T,H) = min Q(H, T;m'). (4.15) 

m' 

The minimization gives: 

j 1 _|_ rri 

zJm + H = -fc^Tln — — — , (4.16) 
2 1 — m 
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Figure 4.3: External-field dependence of the Ising-model magnetization in the 
mean-field approximation. T c is the critical temperature. Solid line: stable states, 
dashed line: metastable states, dotted line: unstable states. 



This is the familiar "equation of state" of magnetic systems. 

The field-dependence of m, Eq. (4.17) is shown in Fig. 4.3. At high tem- 
peratures (T > T c ), the order parameter m behaves as in a paramagnet, this is 
the paramagnetic phase. The spins are partially ordered only because of H. For 
H = 0, the spins are disordered, therefore this phase is also called the disordered 
phase. At low temperature (T < T c ), the spins are spontaneously ordered even 
in the absence of H. This is the spontaneously ordered phase. The system has 
long-range order for J > 0, all the spins are oriented predominantly in the same 
direction, this is the ferromagnetic phase. 

The Gibbs free energy is: 



or: 



m = tanh [f3(zJm + H)\. 



(4.17) 



G(T, H) 
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where m is related to H via the Eq. (4.17). The equation of state is 

(dG(T, H)\ 



M = 



dH 



(M = Nm). 



(4.19) 



The Helmholtz free energy is F(T, M) = G(T, H) + EM: 



F MF (T, M) = Nk B T 



1 + m 1 + m 
In 



1 — m , 1 — m 
H In 



1 



NzJm 2 . 



(4.20) 



The field dependence of the equilibrium Gibbs free energy G is shown in Fig. 
4.4. The equation of state now reads 



H 



fdF(T, MY 
\ dM 



(4.21) 



The Helmholtz free energy F(T, M, N) is shown in Fig. 4.5. 

As before, the critical temperature, where the spontaneous ordering vanishes 
in the absence of external magnetic field H as T is increased, is determined by the 
condition f3 c zJ = 1, 



k B T c = z\J\ 



(4.22) 



In the past, many improvements to the simple MF theory presented here have 
been proposed. With these improvements, the MF critical temperature came closer 
to the exact critical temperature. None of these theories, however, was able to 
bring an improvements to the critical exponents, they kept their mean-field values. 
The real break-through brought the renormalization-group theory, which will be 
introduced in Section 4.3. 



4.1.3 Landau Theory of Continuous Phase Transitions 

In the previous Section we have seen that the trial Gibbs free energy Q(T, H; M') 
was an analytic function which had approximately a parabolic shape as a function 
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Figure 4.4: Equilibrium Gibbs free energy. At low temperatures, when a first- 
order transition at H = takes place, part of the curve corresponds to metastable 
states and part to unstable states. 

of the variational parameter M' above the critical temperature and a double min- 
imum below Tq. Close to the critical point and for H = 0, the minimum of Q is 
close to ml — and we can expand the trial mean-field free energy Q to a very 
good approximation in a power series of m. 



At T > T c , r must be > and below T c , it must be < 0, thus r changes sign and 
in most of the cases it can be approximated by a linear temperature dependence, 



The Landau theory is phenomenological and deals only with macroscopic 
quantities. The theory does not deal correctly with the fluctuations and is a mean- 
field theory. In fact, the Landau theory of continuous phase transitions is much 
more than just simple expansion of Q. It is based on the symmetry of the system. 
Landau postulated that the free energy can be expanded as a power series in m 
where only those terms that are compatible with the symmetry of the system are 
allowed. 




(4.23) 



r = r'(T-T c )/T c . 
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Figure 4.5: Equilibrium Helmholtz free energy. At low temperatures, where the 
transition is first order, F has two minima, corresponding to ±m. The region 
between the two minima is the coexistence region when the system can be partially 
in the +m, and partially in the — m state (domain structure!). 
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Formulation of the Landau Theory 

Critical anomalies take place on long lenthscales. Therefore it is enough and in 
order to change from discrete to continuous description: 

rrii — * m{f) 
hi ->• h(r). 

In formulating the Landau theory let us consider a macroscopic system and 
assume that there exists a (generalized; not necessarily magnetic) order parameter 
density m(r) averaged over microscopic, atomic distances a. Because of the aver- 
age, m(r) is now a continuous function of r (continuous in m and in r)\ There can 
be no fluctuations in m(f) on the lengthscale smaller than a (a^ lattice spacing). 
Correspondingly, m(r) has no Fourier components with wavevector greater than 
k = 2n/a. 

In lattice models, the canonical partition function appropriate to the Gibbs free 
energy of a spin system is: 

Y (T, H) = (4.24) 

{Si] 

(we have to sum over all spin configurations). Now, since m(r) is a continuous 
function of the coordinate r, we must replace the sum by a functional integral over 
all functions m(r). When doing this replacement one must not forget that m(r) 
is equal to the average of a microscopic spin Si over a small volume. That means 
that there can be several microscopic spin configurations which average into the 
same m(r). Let there be Q(m(r)) spin configurations that lead to the same m(r). 
Then, the sum over all spin configurations is equal to the functional integral over 
all functions m(r) where each function m(r) is weighted by Q(m(r)): 

]T -> f Dm(r) Q(m(r)). (4.25) 

The partition function is: 

Y(T,H) = J Dm(r)0(m(r)) e -^w(H;m(r))_ (4 26) 

Using the Boltzmann formula S = ks In Q, we can write Y in the form: 

Y(T,H) = jDmtfe-W&MryTS) 
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= /" Dm(r) e-^ T>JT;m < r » (4.27) 
where we have introduced the Landau free energy: 



Q(T, H; m(r)) = H(H; m(r) - TS(m(r)). 



(4.28) 



Notice the analogy between this equation and the trial Gibbs free energy ob- 
tained previously with the Bragg-Williams method! The Landau free energy is in 
fact a generalization of the trial Gibbs free energy. 

Notice also that the entropy S(m(r)) is not equal to the usual microcanonical 
entropy S (related to the number of configurations in an energy interval AE) 
because we did not make any summation over the energies! S(m(r)) is related 
to the number of microscopic configurations leading to the same function m{r). 
Only after integration over all the states in a certain energy interval we would 
get the usual thermodynamic entropy S. Therefore we can call S(m(r)) a partial 
entropy. 

We still have to work out the functional integral in order to get the (equilib- 
rium) partition function Y(T, H) ! From Y(T, H) the Gibbs free energy is: 

G(T, H) = -k B T In Y(T, H) = -k B T In J Dm e -^(T,ff;m)_ (4.29) 

In the Landau approach we write down the Landau free energy in a most gen- 
eral form as an expansion in terms of the order parameter and its derivatives by 
taking into account all the symmetry-allowed terms up to a given order: 

Q(T, H,m) — J dV (^-rm 2 + um 4 + ■ ■ ■ — Hm 

+ -#|Vm| 2 + ---) . (4.30) 

The highest-order coefficient in the expansion must always be positive, else the 
system were unstable. For the Ising model, this form can also be obtained by 
expanding G, Eq. (4.14), in powers of m. The last term appears if m is site- 
dependent: 



SiSj — > m(r) m(r + a) = m(r) m(r) + m(r) aVm(r) H — a 2 V 2 m(r) + 

2 



(4.31) 
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The first term is already included in the first line of (4.30), the second term van- 
ishes in systems with center of inversion (symmetry argument!) and the last term, 
after partial integration over the whole volume, gives the term proportional to 
| Vm| 2 . The prefactor g must also be positive always. A negative g would lead to 
completely unphysical wild fluctuations in m on extremely short lengthscales. 

Here we have assumed a one-component order parameter. In general, the order 
parameter can have several components or be complex. The coefficients are in 
principle temperature-dependent and depend also on the symmetry of the system. 
For example, the coefficient of the cubic term vanishes if the system is invariant 
under the reversal of m(r) when H = 0. Usually we assume that u is temperature 
independent whereas the quadratic term has a linear temperature dependence, 

r = r't, t = {T-T c )/T c (4.32) 

In the Landau theory we assume that the fluctuations of the order parameter 
are small, that the important values of m lie in a narrow range near its equilibrium 
value that minimizes the Landau free energy. In fact, the maximal contribution 
to the partition function comes from the term with minimal Landau free energy. 
In the thermodynamic limit, the probability of any other configuration will vanish 
rapidly. Therefore we make use of the saddle-point approximation, the functional 
integral (4.29) is equal to the maximum value of the integrand, multiplied by a 
constant. For constant field H, the functional form of m(r) that maximizes the 
integrand (minimizes the Landau free energy), is a constant, determined by the 
conditions 

dg d 2 g 

am om z 

In the saddle-point approximation we immediately obtain G(T, H) from Eq. (4.29): 

G(T, H) = G + min m g(T, H; m) = (4.34) 

In equilibrium, the system is at the bottom of the Landau free energy g. This 
determines m and G(T, H). From G(T, H) we get other thermodynamic quanti- 
ties, as before. 

The advantage of the Landau theory is that it allows one to study a macro- 
scopic system without knowing all its microscopic details. It is based on the sym- 
metry properties of the macroscopic system alone. 
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4.1.4 Critical Exponents in the Landau Theory 

We will now calculate the critical exponents of the Landau model. To be more 
specific, let us investigate the Landau model of magnetic phase transitions, when 

g = G + Qrm' 2 + urn' 4 - m'Hj V (4.35) 

In equilibrium, when dQ /dm! = 0, we get the equation of state: 

rm + 4um 3 - H = (4.36) 

• Order-parameter exponent f3 ( H — 0). For H = 0, and r > 0, the equation 
of state has a single solution, m — 0, there is no spontaneous magnetization 
above T c , as expected. For t < 0, the equation of state has three roots when 
H = 0. The root at m = belongs to an unstable state, so we are left with 
two real roots. The result is: 



r 



m = ±J^-t 1/2 j3 =\ t<0 



k m = £ > 

• Critical isotherm exponent 5 (T — T c ). Now, r = 0, and we have 4wm 3 = 
H, from which it follows 

m=\— =>- 5 = 3 £ = (4.38) 

V 4-u 

• Susceptibility exponent 7. The isothermal susceptibility is defined as 

(dm\ 

From the equation of state, Eq. (4.36), we get: 

r (wj T +12 ""* {m) t - 1 = (4 - 40) 

from which we find: 

— oc r 1 =^ 7 = 1 £ > 

f r'£ 

XT = I , (4.41) 

oc |£| _1 =>-7 = l £<0 
-2r'£ 1 1 
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Specific-heat exponent a (H — 0). The specific heat at constant (zero) field 
is 

<*-t(£) b (4.42, 

For t > 0, m = 0, and G = G = const, (that means that G does not 
depend on the reduced temperature t), so ch is zero. For t < (and H = 0), 

m 2 = — r/Au and 



„/2 

16m 



G = G - —t 2 (4.43) 



It follows: 



=>. a = t > 



c ff = r /2 T c „ „ (4.44) 



8m 



a = t < 



Comment 1: We see that the Landau theory gives mean-field (also called "classi- 
cal") critical exponents. 

Comment 2: At H ^ 0, there is no phase transition, all the quantities vary 
smoothly, only the exponent 5 is defined. 



4.1.5 Short-Range Order and Fluctuations 

Although the correlations and fluctuations are neglected in the MFA, we can still 
get an estimate of the short-range correlations. To obtain the mean-field values of 
the exponents v and r\, which describe the behaviour of the correlation function, 
we have to include the gradient term in the Landau free energy. In the frame- 
work of the mean-field theory, this leads to the Ornstein-Zernike approximation. 
To investigate the correlations, we consider the response of m(r) to a perturbing 
field acting, say, at the origin, H(r) = \8{r) via the "response function", the 
susceptibility: 

m(f) = J d d r'x(f-f)H(f). (4.45) 
The Fourier transform of the last convolution is simply: 

m(q) = x(0)H(q). (4.46) 



86 



I. Vilfan 



Statistical Mechanics 



Chapter 4 



To the lowest order in m and Vm, the Landau free energy is 
g(H;m(r)) = G + J d 



-g\Vm{r)\ 2 + -rm 2 (r) 
+um 4 (r) — Xm(r)S(r) 



(4.47) 



We first make Fourier transforms to the d— dimensional momentum space, 

m(q) = J d d r e _< * r m(f) 



m(r) 



(4.48) 



The Fourier transform of Vm is: 



Vm(f) = 



(2tt) c 



d d q (ye^)m{q) 
^— d J d d q iqe^ r m(q) 



(27T) 

The Fourier transform of a 5— function is a constant 

In the momentum space, the Landau free energy is: 
g(H;m(q))=G + J d d q 



(4.49) 



(4.50) 



1 1 

-gq 2 m 2 (q) + -rm 2 (q) 

+um 4 (q) - Xm(q) 



(4.51) 



Here, we approximated the Fourier transform of m 4 (f) by m 4 (g) ! 

In the mean-field approximation the modes with different momenta q are as- 
sumed to be independent of each other. Therefore we can minimize Q for each 
individual mode separately, 



dQ 



0, 



dm(q) 

and the equation of state in momentum space is: 

(gq 2 + r)m(q) + Aum 3 (q) = A. 



(4.52) 



(4.53) 
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For t > 0, when m vanishes in the absence of the field, we can neglect the cubic 
term, and we get: 

m(q) = —2^—, X(q) = (4-54) 

gqz _|_ r ggz _|_ r 

After a (somewhat tedious) Fourier transform back to real space we find in the 
limit of large r : 

x(r) oc m(r) oc r 2 ~ d exp(— r/£) =>- r\ = (4.55) 

where 

f = \/? ^ = i (4-56) 

Below 7b, m ^ 0, we replace m(f) by m+5m(r), where m is the spontaneous 
magnetization and 8m the response to the infinitesimal perturbing field A. We 
expand in powers of 5m and retain only linear-order terms. We find: 

Sm(q) = / - (4.57) 

and the same exponents as for t > 0. 

In the Ornstein-Zernike approximation (random-phase approximation), the 
modes with different momenta are considered to be independent of one another. 
The free energy is a sum (integral) of contributions from non-interacting normal 
modes. 



4.1.6 Validity of The Mean-field Theory 

The mean-field theories neglect the fluctuations of the neighbouring spins, each 
spin interacts only with the mean value of its neighbours. So, the mean-filed 
theories can only be valid when fluctuations are irrelevant. We will now estimate 
the relevance of fluctuations close to the critical point. We shall calculate the free 
energy of the fluctuations and compare it with the free energy calculated in the 
MFA. 

To the order of magnitude, only those fluctuations can be thermally excited 
whose free energy is of the order ksT (or less). The free energy of a typical 
fluctuation will thus be proportional to k B T. Since the size of the fluctuation is of 
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the order of the correlation length £ (this is how the correlation length is defined), 
the free energy of a typical fluctuation per unit volume is proportional to 

G fluct ex C d oc \t\ vd . (4.58) 

On the other hand, the singular part of the mean-field free energy is governed by 
the specific-heat exponent a: 

G MF oc \t\ 2 ~ a (4.59) 

(we obtain this expression after two integrations of the specific heat over temper- 
ature). For the mean-field theory to be "exact", Gf luct <C G M f in the limit as 

t -> 0, 

|t|-o G M f l*l-»o 

2 — a 

^ud-2 + a>0^d> . (4.60) 

v 

For the Ising model (scalar order parameter) in the mean-field approximation a = 
and v — 1/2 and we get the condition 

d > 4. (4.61) 

The mean-field exponents are correct and the fluctuations are irrelevant only for 
d > 4. d = 4 is called the upper critical dimension of the Ising model. 
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4.2 Monte Carlo Simulations 

With the modern computers, different computer simulations became extremely 
powerful techniques. One sets up in the computer the Hamiltonian describing a 
physical system and tries to see how it behaves. However, an efficient and reliable 
computer simulation is full of tricks and a lot of experience is needed to produce 
reliable results. 

Among different computer simulations, the Monte Carlo (MC) simulations are 
probably the most widely used in statistical mechanics. MC simulations are used 
to study classical as well as quantum systems, continuous spin systems, liquids 
and other disordered systems, gauge theories, etc. Here we shall explain the prin- 
ciples of the MC methods on the prototype example, the Ising model on a square 
lattice. 



4.2.1 Importance Sampling 

In statistical mechanics we calculate the expectation values of thermodynamic 
observables. For a spin system, the expectation value is defined as 



For an Ising system, the sum is over 2 N configurations. This number is extremely 
large for all but very small systems. Therefore one must find a method which 
estimates thermodynamic properties by sampling a small subset of representa- 
tive configurations. One possible strategy would be to scan randomly the whole 
configuration space (scan randomly over all spin configurations). The expecta- 
tion value of A would then be obtained as a sum over the scanned configurations, 
weighted by the appropriate Boltzmann factors. However, this random sampling 
has a very serious drawback. Because of the Boltzmann factor which brings a 
negligible weight to most of the configurations, very few configurations will con- 
tribute to the expectation value of A and a very unreliable estimate will result. 
This problem occurs because only an extremely restricted part of the configura- 
tion space is occupied with a considerable probability in the thermodynamic limit. 
Therefore it makes sense to restrict the sampling only to these states, this is the 
importance sampling. 

Importance sampling is realized if we generate a Markov chain of configura- 
tions, i.e., a configuration is generated from the preceding configuration. The dis- 
advantage of this Markov chain is that the successive configurations are strongly 
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correlated and therefore very long sequence of configurations is needed to obtain 
reliable averages. The expectation value of A is now equal to: 



(A) = lim 



n max n min n=n min 



(4.63) 



where A n is the average of A over n successive configurations. 



4.2.2 Metropolis Algorithm 

In generating successive configurations, the transition probability must be normal- 
ized, ergodic, and it must obey detailed balance. The Metropolis algorithm fulfills 
the above conditions. A new configuration {Sf} is obtained from the previous 
configuration {Si} by flipping one or more spins. In the Metropolis algorithm, the 
probability that the system changes from {Si} to {Sf} is 



r e -p(E f -Ei) if E > E . 

P({S t } {S,}) = 'f (4.64) 



where E { and Ef are the energies of the previous and final states, respectively. 
With this choice of transition probabilities, the system tends asymptotically to a 
steady state in which the probability of a given configuration is e^i 3E d s }\ 

The practical implementation of the Metropolis algorithm is very simple. First, 
a spin is selected either at random or sequentially. The spin is flipped and the en- 
ergies Ei and Ef axe calculated. Finally, a pseudo-random number is used to 
accept or reject the flip with the probability (4.64). The meaning of this condition 
is straightforward. If the systems gains in energy by flipping the spin, the flip 
is accepted. It is also accepted if it costs energy to flip the spin, but the energy 
difference is small so that the Boltzmann factor is larger than the random number 
(analogy to thermal excitations!). The flow chart of a MC program is shown in 
Fig. 4.6. 



4.2.3 Practical Details and Data Analysis 

• Initial Configuration. In principle, the initial configuration should not in- 
fluence the final configuration. However, the equilibrium may be reached 
faster if a proper initial configuration is chosen. An example, showing the 
approach to thermal equilibrium in a d = 3 Ising ferromagnet, is shown in 
Fig. 4.7. 
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Set up the lattice 
(Initial configuration) 


1 


f 



Choose a site 
Flip the spin 
Calculate b = exp(-AE/kT) 
Choose a random number 
< r < 1 



b> 



Accept the flip 



k = k + 1 



b <r 



Reject the flip 



Calculate the physical quantity 
of the configuration k 



Calculate the average of 
the quantity A: 
A — 1 V* 2 A, 



Figure 4.6: Flow diagram of a Monte-Carlo simulation with Metropolis algorithm. 
During the first k\ steps the system still approaches to its equilibrium configura- 
tion, therefore these configurations are not included in calculating the average 
properties. 
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0.2- 



30 000 60 000 90 000 120 000 150 000 

Number of configurations 

Figure 4.7: Relaxation of magnetization towards its equilibrium value in a Monte- 
Carlo simulation of an Ising system on a 10 x 10 x 10 lattice with periodic boundary 
conditions. Solid lines are "instantaneous" magnetizations and dashed lines their 
equilibrium estimates. (From: K. Binder and Z. Rauch, Z. Phys. 219, 201 (1969).) 



93 



I. Vilfan 



Statistical Mechanics 



Chapter 4 



3.0 



1.0 



0.3 



0. 1 



T<T t 

slope.! 



» N ■ 

oN : 

• N : 

- <N : 
□ N i 
a N : 
+ N i 

• N ' 

i 

0.1 



10 

1 4 
20 
24 
30 
40 
50 
60 



•An 



T>T C 



slope -| 



1.0 
x = N 1/v t 



\ 



10.0 



Figure 4.8: Scaled magnetization vs. scaled temperature of an Ising model on a 
square lattice. (From: D.P. Landau, Phys. Rev. B 13, 2997 (1976).) 

• Boundary Conditions. Different boundary conditions are used. In periodic 
boundary conditions the system is periodically repeated in all directions to 
mimic an infinite system. Alternatively, one may describe the effect of the 
surrounding spins on an (inevitably) finite lattice by a mean field, acting on 
the spins at the border. Open boundary conditions, when the finite lattice is 
isolated and the spins at the surface have less neighbours than other spins, 
have also been used. 

• Finite-size effects. Lattices, considered in MC simulations, must have a 
finite size (L) whereas we are usually interested in the macroscopic prop- 
erties. The question is, how to take the limit of L to oo. Usually, one 
calculates an average for several different lattice sizes and then extrapolates 
the average to infinite L. This extrapolation is simple away from the critical 
point where the correlation length is small compared to L. The problem 
becomes non-trivial in the critical region because £ which would normally 
diverge, cannot exceed L. As a consequence, singularities associated with 
the critical point are rounded. If possible, scaling forms should be used in 
extrapolating L — >• oo (see Fig. 4.8). 

• Random-number generator. Generation of uncorrelated random numbers is 
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not a trivial problem. Computer languages have built-in "pseudo-random- 
number" generators. Not always these pseudo-random-number generators 
are random. In particular, if random sampling is used, then it might happen 
that some points on the lattice are visited more often than the other (which 

- in the worst case - may never be visited). Pseudo-random-number gen- 
erators use algorithms to calculate a sequence of (random) numbers. Care 
must be taken to initialize random number generators with different "seeds" 
in order not to obtain the same sequence of random numbers each time. 

• Statistical errors. To obtain reliable results, an average must be taken over n 
much larger than the number over which the MC states are correlated. This 
becomes more difficult near the critical point because of critical slowing 
down. 

• Free energy. Monte-Carlo simulation samples only a very small part of the 
phase space whereas to obtain the partition function or the free energy, we 
must sum over the whole configuration space. In MC simulations, there- 
fore we calculate the free energy by integrating one of its derivatives, e.g., 

— / MdH. This means that we need to calculate M(H) in a range of H . 

• Conservation laws. Above, the spins were allowed to flip uncorrelated from 
one another. This works fine for magnetic systems, where the order param- 
eter is not conserved. For binary alloys, e.g., however, the number of atoms 
is conserved and so is the order parameter. So, opposite spins have to flip 
simultaneously (one "up"^"down" and another - usually a neighbouring - 
one "down"^"up"). The conservation laws have important consequences 
for the kinetics of fluctuations (slowing down). 
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4.3 The Renormalization Group Methods 

Some typical results of MC simulations of a d = 2 Ising ferromagnet are shown 
in Fig. 4.9. At high T (T = 2T C ), there is only short-range order, the spins 
form small clusters. The correlation length (approximately equal to the linear size 
of the largest cluster) is small. Close (but above) T c , somewhat larger patches 
in which most of the spins are lined up in the same direction begin to develop. 
When the system reaches T c , these patches expand to infinite size, but fluctuations 
of smaller scale persist! At the critical temperature, spins form clusters at all 
lengthscales, including one infinite-size cluster (which cannot be seen on a finite 
system, but we know that the correlation length has to diverge at T c ). As a result, 
all scales of length must be included in a theoretical description. Another very 
important conclusion is reached if we substitute a block of 3 x 3 spins by a single 
spin with the orientation prescribed by the majority rule (block spin "up" if the 
majority of the spins in the block is "up", and vice versa). In this way the length 
scale of the lattice is changed by a factor 6 = 3 each time. This is a "real-space 
block-spin renormalization-group transformation". 

Examples of such transformations are shown in Figs. 4.10-4.12. Below Tc, 
the system has long-range order but, because the system is close to T c , the fluc- 
tuations form large clusters of oppositely ordered spins. However, under each 
transformation the size of the unit cell increases, the fluctuations become smaller 
and smaller and those fluctuations, which are smaller in size than the lattice spac- 
ing, disappear. Under successive transformations, thus, the system becomes more 
and more ordered, it effectively approaches zero temperature, see Fig. 4.13. 

Above T c (Fig. 4.11), there is no long-range order, spins form random "up" 
and "down" clusters. Under each transformation, the correlation length (expressed 
in units of the new unit cell) decreases, the clusters become smaller and smaller 
as if the temperature were higher, but they never disappear. At T = Tc (Fig 
4.12), the clusters are of all length scales, including one infinitely large cluster. 
The correlation length is infinite (exceeds the size of the lattice used in the MC 
simulation). Under successive transformations, it remains "infinite", there is no 
change in the distribution of cluster sizes, the system is invariant under scale trans- 
formations ! The system remains at the critical temperature and the critical point 
is a. fixed point. Other (trivial) fixed points are T = and T = oo because the 
system away from T c tends to one of these points under successive scale transfor- 
mations. Decrease of the correlation length under successive transformations is an 
extremely useful property because it makes the fluctuations uncorrelated and we 
can solve the system using an approximate (e.g., mean-field) theory to calculate 
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Figure 4.9: Results of MC simulations of a d = 2 Ising ferromagnet at different T. 
At high T (T = 2T C ; top), the correlation length is small. Close (but above) T c 
(middle), somewhat larger clusters begin to develop. When the system reaches T c 
(bottom), these clusters expand to infinite size, but fluctuations on smaller scale 
persist! (From: K.G. Wilson, Sci. Am. 241, 140 (1979).)) 
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Figure 4.10: Block-spin transformation is applied to a lattice repeatedly at T = 
0.99T C . Each time the number of spins is reduced by a factor of 9, elucidating 
the behaviour of the system at a larger scale. The correlation length decreases 
under successive transformations. Clearly seen is the decrease of the correlation 
length and a suppression of the fluctuations. Under each transformation the system 
appears more ordered, the system flows under the renormalization transformations 
towards zero temperature. (After K.G. Wilson, Sci. Am. 241, 140 (1979)) 
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Figure 4.11: As Fig. 5.9 but for T = 1.22T C . Under each transformation, the 
correlation length decreases, the system flows towards infinite temperature. (After 
K.G. Wilson, Sci. Am. 241, 140 (1979)) 
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Figure 4.12: As Fig. 5.9 but for T = T c : The correlation length is infinite 
and remains infinite under successive transformations, the system remains at the 
critical temperature. (After K.G. Wilson, Sci. Am. 241, 140 (1979)) 
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Figure 4.13: Block-spin renormalization at H = 0. Upon successive decimations 
of an initial configuration at T < T c the system becomes more and more ordered 
as if the temperature were lower and lower. Eventually, the system becomes com- 
pletely ordered, this is the T = fixed point. Upon decimations of a state at 
T > Tc, the system becomes more and more disordered, the correlation length 
decreases as if the temperature were higher, until the correlation length becomes 
on the order of the lattice spacing. This is the infinite-temperature fixed point. 
A system at T c stays, upon successive decimations, at T c . This is another fixed 
point, the critical point. 

the properties in that region and then transform back to the original lattice. This 
is the basic idea of the renormalization-group methods and we shall exploit this 
idea on some simple examples. 

4.3.1 A Trivial Example: The d = 1 Ising Model 

In the renormalization-group methods the temperature changes under successive 
transformations, therefore it is convenient to work with the reduced Hamiltonian, 
we divide H by ksT. The reduced Hamiltonian of an Ising Hamiltonian with NN 
exchange interaction in an external field is 



H = -K y £S i S i+1 -h , £S i , K=— h=—. (4.65) 
a a k B T k B T 



Recursion relations. If we want to make several successive renormalization- 
group (RG) transformations, the structure of TC must not change after each trans- 
formation. We are thus seeking for such a transformation from the old Hamil- 
tonian TC to the new Hamiltonian TC' that keeps the structure of the Hamiltonian 
unchanged under successive transformations. In addition, the transformation must 
be such that it does not change the partition function. That is necessary because 
this ensures that the physical properties of the system are not changed. The parti- 
tion function of the one-dimensional Ising model in an external field is written in 



J 



H 
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Figure 4.14: Decimation on the d — 1 Ising model. 

the form 

Z(N,K,h) = Y,e- H = J2e KEiSiSi+1+hEiSi 

{S} {S} 
_ ^2 e K(SiS 2 +S 2 S a )+hS2 + ^(Si+S 3 ) 

{S} 

y^QK(S-iS4+S4S 5 )+hS4+^(S3+S5) x . . . (4 

We divide the lattice into the sublattices x and o, denote the spins on the x sub- 
lattice by <Ti (Fig. 4.14), 

Z = ^^ e <T2[K(S 1 +S 3 )+h]+%(S 1 +S 3 ) e *4[K(S 3 +S s )+h]^(S 3 +S B ) 

{S} M 

(4.67) 

and carry out a partial trace, i.e., we sum over 2 Ar / 2 configurations of the {a} spins 
sitting on the x sublattice: 

Z(N, K,h) = [e^+tX^m + e (-tf+!)(Si+s s )-h] x . . . . (468) 

{S} 

The entire partition function must not change and H must keep its structure, there- 
fore we have the condition 

Z{N,K,h) = e N ^ K ' h) Z{-,K',h') =c N9 Y,^ H ' ( 4 - 69 ) 

2 {s} 

where 

V! = -K' £ SiS i+2 - ti Si (4.70) 

jgodd ieodd 

and g is a constant. This will ensure that the physical properties will not change 
under successive iterations. The two conditions (4.70) and (4.69) imply that 

e (K+^)(Si+S i+2 )+h + Q (-K+\)(S l +S l+2 )-h 
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= e .K'S l S I+ 2 + ^(S i +S i+2 )+2g _ (4.71) 

The last equation has to hold for any S i: S i+2 — ±1: 

e 2JsT-2h + e -2X = Q K'-h'+2g 

& +h + e -h = Q ~K'+2g^ (4?2) 

from which we get: 

^2K+2h + ^ e 2A"-2ft + Q -2K^ 

— 



(e A + e- ft ) 2 



e 2fe ' = 



e 2in-2/t _|_ e -2X 
& 2K-2h _|_ £ -2X 



e 8 fl = ^ e 2K + 2h + e -2Ar^ e 2Ar-2/, + e -2^^ + e -^ 2 (473) 



or: 



1 cosh(2X + /i) cosh(2fT - /i) 
i\ = — m 



4 cosh /i 

, , , 1, cosh(21f + h) 

h = h+ 2 ln cosH2K-h) 

g = I In [16 cosh(2X + h) cosh(2iT - h) cosh 2 /il . (4.74) 
8 l J 

Equations (4.73) or (4.74) are the recursion relations and determine the fixed 

points and the flow diagram of the system. 

In each iteration, the number of degrees of freedom is reduced by one half, the 

new Hamiltonian H has only one half of the previous spins, 

N' = N/b 6 = 2, (4.75) 

and the lattice spacing is increased, 

a = ba. (4.76) 

Other lengths which are measured in units of the lattice spacing are reduced by a 
factor b. For example, the new correlation length is 

e = i/b. (4.77) 

The remaining spins on the decimated lattice interact with their new nearest neigh- 
bours through the renormalized coupling constants K' and are subject to renor- 
malized fields h! . 
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Fixed points. At a fixed point the parameters (the coupling constants - in our 
case K and h) do not change under successive decimations. This means that the 
system at a fixed point stays at this point in the parameter space. With the new 
variables 

x = e" 4 * y = e' 2h z = e"^ 

0<(x,y,z)<l (4.78) 

the recursion relations become 



X = X 



y' = y 



(x + y)(l + xy) 
x + y 
1 + xy 



z ' = z 2 xy 2 - 1 —— - (4.79) 

(x + y){l + xy){l + y) 1 

The recursion relations for x and y do not depend on z. Physically, this means 
that the singular behaviour of the free energy does not depend on a shift in the 
energy scale. Therefore, we first investigate the flow diagram and fixed points 
in the (x, y) plane, see Fig. 4.15. From the first equation it is easy to see that 
x' — x if x — 1, i.e., x = l is a fixed point of this system for arbitrary y . We 
denote the fixed points by a *. There is a line of fixed points at x* — 1 for 
< y* < 1, corresponding to an infinite-temperature sink. Along this line, the 
spins are completely disordered and £ — > 0. If we set x — 0, the second equation 
reduces to y' = y 2 , which has two fixed points, one at y* — 1 and one at y* = 0. 
The fixed point at (x*, y*) = (0, 1) (this means H = and T = 0) corresponds to 
the critical point where £ — > oo. This is correct because we know that the d — 1 
Ising model has Tc = 0. This fixed point is unstable and any point in its vicinity 
will flow under RG transformations away from it towards the fixed points on the 
line at x* = 1. The critical behaviour of the d = 1 Ising model is governed by 
the least stable fixed point at (x*,y*) = (0, 1). The point at (x*, y*) = (0, 0) is 
another fixed point which corresponds to T = and H = oo. At this fixed point, 
the spins are completely ordered. 



4.3.2 RG Transformations on the d = 2 Ising Model 

The recursion relations of the d — 1 Ising model could be formulated exactly. 
Usually, however, farther-range interactions are introduced under successive dec- 
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Figure 4.15: Schematic flow diagram of the d — 1 Ising model. (After J.M. 
Yeomans, Statistical Mechanics of Phase Transitions, Oxford 1992.) 



imations and approximations have to be made in order to get the closed-form 
recursion relations. 

As an example, let us apply the decimation on the d = 2 Ising model on 
a square lattice in the absence of an external field. We start with NN (nearest 
neighbours) exchange interactions and sum the partition function over all possible 
configurations of the a spins (these are the spins that will be decimated, they are 
denoted by x in Fig. 4.16). 



{s} {s} m 

_ J~J ^QK(Si+x+Si- x +S i+ y+Si-y) _|_ ^-K(Si +x +Si~ x +Si+y+Si- y ) 

{S} i 



(4.80) 



S± x and S± y are the "right/left" and "upper/lower" neighbours of <jj, respectively. 
The conditions are that under decimation Z must not change and that H should 
have the same structure as the original Hamiltonian. However, we soon realize 
that we cannot fulfill the above conditions with the renormalized form of the initial 
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Figure 4.16: Decimation on the d = 2 square lattice. The spins denoted by x 
are decimated. K is the NN interaction and L the NNN interaction. After deci- 
mation, the lattice constant increases by b = and previous NNN become NN 
(nearest neighbours). Solid lines represent the original lattice and dashed lines the 
decimated lattice. 

Hamiltonian. After decimation we get a more complicated Hamiltonian: 

7~L = — ^ [K (Si + X Si + y + Si + ySi^ X + Si^ X Si-y + S^y (S'j+j; ) 

i 

H~ L (Si+ X Si — X + Si+ySi — y) + Q S i+Q^Si+yS i — X S i — y + C ] , 

(4.81) 

where X' is the new NN interaction, L' the new NNN (next nearest neighbours) 
interaction, Q' the four- spin interaction, and C a constant. Only with the Hamil- 
tonian of this form, it is possible to satisfy the condition Z' = Z for arbitrary 
configurations of S± XtV . Although we started with only NN interactions, decima- 
tion generated NNN interaction and even four-spin interactions. So, we have to 
start from the beginning again with a more complicated Hamiltonian which in- 
cludes all these interactions. In the following we shall neglect Q' and assume that 
L is small. In this case we obtain the following approximate recursion relations: 

K' « L + 2K 2 

U « K 2 (4.82) 
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Olx^O 



Figure 4.17: The renormalized interaction after a decimation is L + 2K 2 . 



In fact, these relations could be guessed immediately, see Fig. 4. 17. The new near- 
est neighbours are the previous NNN, therefore the term L in the first equation. 
Besides, they were connected by two pairs of consecutive bonds, this brings IK 2 . 
In the second equation, the new next nearest neighbours (NNN) are connected by 
two consecutive NN bonds on the old lattice - each of them contributing a factor 
K. 

The recursion relations 4.82 have two trivial fixed points, (K*, L*) = (0, 0) 
and (K*,L*) = (oo, oo) and one non-trivial fixed point at (1/3, 1/9). The first 
one is the infinite-temperature, paramagnetic fixed point, the second one is the 
zero-temperature, ferromagnetic fixed point. The non-trivial fixed point corre- 
sponds to the critical point of the system. 



4.3.3 General Properties of the RG Transformations 

So far, we have seen that the renormalization-group transformations work by 
changing the length scale and by reducing the number of the degrees of freedom. 
Only at criticality the system remains unaltered under the change in scale and the 
critical behaviour is determined by a fixed point of the transformation. 

We start with an effective Hamiltonian, which is characterized by a set of 
coupling constants and fields K, L, ■ ■ ■. Each set of these constants defines a point 
in the "Hamiltonian space," the space of coupling constants. The renormalization- 
group transformations work in this space, 



w = nn. 



(4.83) 



The renormalization-group operator 1Z reduces the number of degrees of freedom 
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from N to N'. This defines the scale factor b: 

N' 



(4.84) 



In the previous examples, the reduction was made in real space either by block- 
spin transformation or by decimation; later we shall do the reduction also by inte- 
grating out large momenta (short wavelengths) in momentum space. 

The essential condition to be satisfied by any RG transformation is that the 
partition function must not change, 



Zn'(7~L') = Z N {H). 



(4.85) 



Therefore the total free energy does not change but the free energy per unit cell 
(spin) increases as 



W) = b d f(H). 



(4.86) 



All lengths, measured in units of the new lattice spacing, are reduced by the factor 
b. Thus, the correlation length scales as 



(4.87) 



Fixed points. The iteration of 1Z traces a trajectory in the Hamiltonian space, 
this is called the "Hamiltonian flow." The trajectories end in fixed points H*, de- 
fined by 

h* = nn* n' = n = n*. (4.88) 

At the fixed points, the system is invariant under subsequent scale changes and 
thus the correlation length does not change 

= £ = (4.89) 

Let us now take a closer look at the Hamiltonian flow near the fixed points. 
We introduce a reduced Hamiltonian which we write in the general form 

H = pi-$ (4.90) 

and investigate the flow in the parameter space spanned by the generalized fields 
jl, conjugate to the operator <fi which contains product of operators (like spins). 
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• Example: In the case of the d = 2 Ising model on a square lattice, the 
components of include the terms linear in S (the conjugate filed is h), 
quadratic in S (coming from interactions with NN, NNN neighbours), four- 
spin terms, etc. The corresponding conjugate fields are the magnetic field h 
and the coupling constants K, L, etc. In the beginning, many components 
of jl are zero, but under successive transformations, more and more complex 
interactions emerge and an increasing number of components of jl becomes 

7^0. 

Under successive RG transformations, the Hamiltonian evolves and the system 
moves through the parameter space: 

jl' = Kj2 (4.91) 

until it reaches a fixed point when: 

/T = /r=/r. (4.92) 

Every particular type of criticality corresponds to a particular fixed point and de- 
viations from the critical point are related to deviations from the corresponding 
fixed point. To get insight into the critical properties we must investigate the flow 
in the Hamiltonian space in the vicinity of the fixed points. 
Near a fixed point we linearize the operator 1Z. We write 

jl = jT + 5j2 

jH = jT + SjH (4.93) 
where the small deviations from jl* are related by 

on 



5jl = jl* + C(jT)5jl (4.94) 



6 ft = C(jT)5jl. (4.95) 

The linear matrix C is evaluated at the fixed point and is thus a constant. The 
eigenvalues A; and the eigenvectors Vi of C, 

CVi = (4.96) 

determine the critical properties of the Hamiltonian. 
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The eigenvalues A« are functions of b. Two successive transformations with 
scale change b make a total scale change b 2 , they are thus equivalent to a single 
transformation with the scale change b 2 . 

Xi(b)Xi(b) = \ t (b 2 ) (4.97) 

and the eigenvalues must be of the form 

Xi(b) = b x \ (4.98) 

Xi are the critical exponents and, as we shall see, they are closely related to the 
usual critical exponents a, (3, ■ ■ ■. So, if we know A;, we can obtain the critical 
exponents Xi from 

• Example 1: For the d — 1 Ising model, jii = x and \ii = y. After lineariza- 
tion of the recursion relations (4.73)around the fixed point at ft* = (0, 1) we 
get: 

5x' = 45x 5y' = 25y. (4.100) 

Hence, C is diagonal in this representation. The eigenvalues and corre- 
sponding eigenvectors are: 



Ai = 4 v x 



A 2 = 2 v x = ( ° ) . (4.101) 



1 

It follows that the two critical exponents are both positive: 

Xl = 2 x 2 = l. (4.102) 

Example 2: For the d = 2 Ising model, /ii — K and fi 2 = L, and after 
linearization Eqs. (4.82) become 

SK' = AK*5K + 5L 

5L' = 2K*5K. (4.103) 
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For the non-trivial fixed point at (K*, L*) = (|, |), the RG transformation 
matrix is: 



(4.104) 



AK* 1 




4/3 1 


2K* 




2/3 



The eigenvalues of C are: 

\i = 2 -±^ „ 1.72, 
3 

and the eigenvectors: 



Ao = 



1 



10 



2.58 
1 

-0.58 
1 



-0.387 (4.105) 



(4.106) 



(Question: what does negative A 2 mean physically, how does the system 
behave under successive RG transformations if A is negative?) The corre- 
sponding critical exponents are: 



x 1 



_ 2 ln[(2 + ^l0)/3] 



In 1(2 



x 2 



In 2 

10)/3| 



In 2 



57 



-2.74, 



(4.107) 



one is positive and the other one is negative. 



Flow in the parameter space. We first express 5 fx near the fixed point at ft in 
terms of the eigenvectors Vi of C. In case of the d = 1 Ising model this means that 



In general we write 5/2 as: 



5x' — X±Sx — 45x 
by 1 = \ 2 5y = 2Sy. 



(4.108) 
(4.109) 



The coefficients hi are the linear scaling fields. They are the components of the 
generalized field in the coordinate system of the eigenvectors. Under a renormal- 
ization 



5ft = LSfx = ^2 hiCUi = hKvi, 



(4.110) 
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where Aj = b x \ On the other hand: 



5jl' = h'iPi 



(4.111) 



therefore 



K = b Xi hi 



(6>1). 



(4.112) 



It follows that the flow of the Hamiltonian in the parameter space close to the fixed 
point depends crucially on the exponents x^ For positive Xi (A; > 1), the scaling 
filed hi increases under repeated RG transformations and drives the system away 
from the fixed point. In this case, hi is a relevant scaling field. If the exponent 
%i is negative, the corresponding scaling field hi decreases under subsequent RG 
transformations, thus driving the system towards the fixed point. Such fields are 
called irrelevant scaling fields . 

Thus the stability of the fixed point and the topology of the flow lines around 
it depend on the number of relevant and irrelevant scaling fields associated with 
it. To illustrate this, Fig. 4.18 shows the the trajectories in the parameter space 
with two irrelevant and one relevant scaling field. 

In general, any system which has initially the component of the relevant field 
different from zero, will after repeated RG transformations eventually be driven 
away from the fixed point (dotted lines in Fig. 4.18). Physically this corresponds 
to the system moving away from criticality. 

The points in which the relevant field vanishes, define the critical surface. 
All the points on the critical surface flow into the fixed point. The points on the 
critical surface represent systems with an infinite correlation length. Under RG 
transformations, the systems stays on the critical surface (£ remains infinite), but 
it moves towards the fixed point. 

• Example: The d — 1 Ising model has hi — x and hi = y. For the non-trivial 
fixed point at x* = and y* — 1 the critical exponents are x\ = 2 > and 
x 2 = 1 > 0, that means that both, temperature and the external magnetic 
field, are relevant. 

In the case of the d = 2 Ising model on the square lattice, x\ > and 
x 2 < 0, as we have seen. It has one relevant scaling field (h{) and one irrel- 
evant field (h 2 ). The relevant scaling field is related to reduced temperature 
whereas the irrelevant field is related to L, the coupling with NNN. Notice 
that here we considered the Ising model without external field. If the field 
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h. 
































Fixed _ , - 
point * 












' Critical 


\ 


\ 




1 


surface 



Figure 4.18: A fixed point with one relevant (hi) and two irrelevant (h 2 and hs) 
scaling fields. All the points on the two-dimensional critical surface flow to the 
fixed point (solid trajectories). The points slightly above (or below) the critical 
surface first flow towards the fixed point but eventually they are driven away from 
it as the relevant scaling field increases (dotted trajectories). 



113 



I. Vilfan 



Statistical Mechanics 



Chapter 4 




(c) 

Figure 4.19: Topology of the RG flow near fixed points with: (a) two relevant 
fields, (b) two irrelevant fields, and (c) one relevant and one irrelevant field. 

is included, then we have two relevant scaling fields and one irrelevant (the 
magnetic field is always relevant). 

Figure 4.19 shows possible scenarios in the case of two scaling fields. Ferro- 
magnets have two relevant fields, external field and temperature (or NN exchange 
interaction) and the fixed point is of the type shown in Fig. 4.19 (a). 

Universality. The critical behaviour is completely determined by the behaviour 
of the system close to the fixed point, where the transformation matrix can be 
linearized. The transformation matrix depends on the values of the scaling fields 
at the critical point and not on the initial values of the irrelevant scaling fields. 
Thus, all systems which flow close to the fixed point, independent of where in the 
parameter space they start, will exhibit the same critical exponents determined by 
the eigenvalues of the linear transformation matrix at that fixed point. The critical 
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exponents are universal and do not depend on the initial point in the parameter 
space (like the coupling strength). 



Critical exponent v. To illustrate how the critical exponent v is obtained from 
Xi, let us consider ferromagnetic systems. We know that there the critical point is 
specified by H and T. These are the relevant fields for magnetic systems. There- 
fore we identify hi by the reduced temperature, hi — t and hi by the reduced 
magnetic field, hi = h. Under RG transformations £ and t scale as 

g = T t' = b Xt t (4.113) 

b 

The two equations tell us that £ Xt t is constant, it is invariant under RG transfor- 
mations. Therefore 

£(t)cxt- 1/xt (4.114) 
Since, by definition of v, £ oc t~ v , we get 



1 

v — — 
x t 



(4.115) 



The last equation can be obtained also in a different way. Unless the system is 
initially at the critical point, the correlation length will decrease under successive 
block-spin RG transformations. The system will eventually reach a point where 
the correlation length £' will become equal to the size of the renormalized unit cell, 
so that the nearest-neighbour spins become uncorrelated. This point corresponds 
to a certain temperature t which cannot depend on the initial temperature. Sup- 
pose, the system is initially at the temperature t, and arrives at t after n iterations. 
Then, 

t = b nxt t and £ = l = ^/b n ^ i = b n (4.116) 

From these equations we get Eq. 4.1 15 again. 

• We return again to the d = 2 Ising model. Above, we considered the model 
at h — 0. Thus, the exponent x t corresponds to the relevant exponent x\ = 
1.57 from Eq. (4.107) and therefore v = 1/xi = 0.64. This is better than 
the mean- field value of 1/2, but still far from the exact value v — 1. The 
discrepancy has its origin in the approximations we have made to obtain the 
RG recursion relations. 
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Other critical exponents. The original and renormalized Hamiltonians Ti and 
Ti! correspond to positions in the parameter space, jl and jl' . We apply relation 
(4.86) to the singular part of the reduced free energy per spin: 

9s ([l') = b d g s (fl). (4.117) 

Near a fixed point jl and jl' can be expressed in terms of the scaling fields hi and 
h\ which are related by (4. 1 12). Hence eq. (4. 117) becomes: 

g s (h 1 ,h 2 ,...) = b- d g s (b^h 1 ,b X2 h 2 ,...), (4.118) 

where g s is a homogeneous function of jl. This is the scaling form of the singular 
part of the free energy. This equation holds for arbitrary b because we can repeat 
transformations many times. In fact, b can be any real number > 1! Again we 
identify h± — t, h 2 — h and assume that other scaling fields are irrelevant. Then, 
eq. (4.118) becomes 

g s (t,h) = b- d g s (b x H,b^h). (4.119) 
Now, to come to the specific-heat critical exponent a, we remind that 

ch = (j^ = gtt(h = 0) oc (4.120) 

Differentiating Eq. (4. 119) twice with respect to temperature and putting the scal- 
ing field h equal to zero gives: 

dr 

gtt(t,h = 0) oc b- d g TT (r,0)— = r d+2 ^ TT , (4.121) 

where we have used r = b Xl t. Initially, t is small whereas b > 1 so that the 
product b Xl t increases with each iteration. We repeat the iterations until 

b Xl \t\ = l. (4.122) 

In this way, the arguments of g TT and therefore also the function g TT are constant, 
all the temperature dependence of the specific heat is in the prefactor h~ d+2xi : 

g tt (t, h = 0) oc \4 d ' 2x ^g a (±l, 0) (4.123) 

and 

a = 2 - — . (4.124) 

Xi 
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Similarly we find 

Xi 

2x 2 — d 

7 = 



5=^— (4.125) 

d — x 2 

We have expressed five critical exponents in terms of two independent expo- 
nents x\ and x 2 . That means that there must be relations between the exponents 
a - • -v. One can easily verify the scaling laws: 

a + 2f3 + 7 = 2 
7 = (3(5 - 1) 
2-a = dv. (4.126) 



4.4 Momentum-Space Renormalization 
Group 

RG methods in momentum space are much more powerful and more widely used 
than their counterparts in real space. The problem is only that they are much 
more complicated than the above real-space prototype examples. We shall define 
the RG transformation for the Landau model, that means that we shall now use 
continuous instead of discrete spin variables. The main reason is that discrete 
spins become essentially continuous after (block-spin) transformations. 
We start with the free energy of the form: 



H(h,m(r)) 



v 



g\Vm[r)\ 2 — m(r)h 



+ -rm 2 (r) _|_ Mm 4 ( r ) 



(4.127) 



This is also called the Landau-Ginzburg-Wilson (LGW) Hamiltonian. The main 
procedure of all the RG transformations is "coarse graining"(i.e., transformation 
to larger unit cells) followed by rescaling (to make the system look like the original 
one). In real- space RG this was made either by introducing block spins or by 
decimation. Now we shall make coarse graining in momentum space. Instead 
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of expanding the size of the unit cell in real space, we shall shrink the size of 
the Brillouin zone in the momentum space. In the following the method of the 
momentum-space renormalization group will be demonstrated on the Gaussian 
model. 

4.4.1 The Gaussian Model 

The Hamiltonian of the Gaussian model in the momentum space is: 

K= \ Y,(r + 9q 2 )\m(q)\ 2 -hm(0). (4.128) 

kKA 

We assume a homogeneous magnetic field and we are only interested in the stabil- 
ity of the m = solution (r > 0), therefore we take u = 0; this will simplify the 
summation over q). A is the cutoff momentum (Brillouin zone boundary). Taking 
u = means that the modes with different momenta become uncoupled. 

The corresponding partition function is equal to the functional integral over al 
functions m(r). In momentum space, the functional integral becomes: 

Dm(r) -»• TT / dm(q) (4.129) 

and the partition function: 

Y(t,h)= FT f°°dm( g )e-^w<A (r+3<,2)|m(<!)|2 - ta(0) . (4.130) 

kl<A 7 ° 

The RG transformation consists of "integrating out" the wavevectors that are 
outside the sphere of the radius A/b and then rescaling. One RG cycle consists of 
three steps: 

1. Integration. After "integrating out" the momenta within a thin shell be- 
tween the radii A/b and A, the partition function (which should not change) 
becomes: 

Y = e n TT fdm(q)e- n ', (4.131) 

\q\<A/b J 

where the first term, exp(fi), is a function of A, b, and coupling constants. It 
comes from the integration over momenta in the shell between A/b and A. 
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The second term is the partition function of the states with momenta \q\ < 
A/b. Notice that this separation into two contributions was only possible 
because we set u — 0, the partition function is then a product of independent 
Gaussian integrals. The new Hamiltonian, 

H'(m(q))= ]T hr + gq 2 )\m(q)\ 2 -hm(0), (4.132) 

\q\<A/b Z 

depends only on m(q) with \q\ < A/b and has the same structure as the 
original one. 

2. Rescaling. The second step restores the old cutoff by blowing up the radius 
of integration to the original value A in such a way that we change the 
variable of integration to 

q' = bq. (4.133) 

In this way the cutoff momentum is restored back to A. This step corre- 
sponds to putting a' = ba in real-space RG. The Hamiltonian now reads: 

H'{m) = b- d £ I( r + ^)|m(()| 2 -M0). (4-134) 

|<?'|<A 1 

3. Normalization. The step 3 restores the standard normalization of the order 
parameter. This step corresponds to putting the new spin equal to 1 after 
block-spin or decimation in real-space RG. Renormalization is restored with 
g' — g if we set: 

m'(q') = (4.135) 

The transformed Hamiltonian and the partition function (which must not 
change) now look like 

H'(m)= ]T Ur' + gq' 2 )\m'(q')\ 2 -h'm'(0). (4.136) 

k'KA 

and 

Y = e n U I dm(q')e-P H ' . (4.137) 
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In this way we integrated out large momenta and transformed H to its orig- 
inal form. The recursion relations are: 

r' = b 2 r 

h' = b^' 2 h. (4.138) 

Since b xs 1 we can write b = 1 + dr, thus, 

b n « 1 +n\nb= 1 + ndr. (4.139) 

We can say that r counts the number of RG cycles. Then, r and h obey the 
following differential equations: 

dr 
d7 
dh 
frr 

with the solution 



2r 

d + 2 



h 



(4.140) 



r = r e 
/i = h e XhT 



x t = 2 



x h 



d + 2 



(4.141) 



Since q is continuous (for macroscopic systems), the RG transformations was 
carried out in infinitesimal steps (b w 1) and the problem was reduced to solving 
a set of differential equations. 

The Gaussian model has two trivial fixed points at h* = 0, one at r* = oo 
(very high T, m — 0) and one at r* = — oo (very low T, where m — > oo because 
we set w = 0!), and one nontrivial fixed point (the critical point) at h* — r* — 0. 
At the nontrivial fixed point the critical exponents are both positive and r and h are 
both relevant scaling fields. From we shall calculate other critical exponents 
in the Section 4.4.3. 

As we shall see in the next Section, for d < 4 the Gaussian model is unsta- 
ble with respect to the parameter u coming from the four- spin interaction in the 
Hamiltonian. 
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4.4.2 The Landau- Wilson Model 

In the Gaussian model, we have neglected the fourth and higher-order terms in the 
Landau Hamiltonian. The next step is to include also the m 4 term. This will lead 
to non-MFA exponents, as we shall see. Here we shall bring only an outline of 
the RG procedure and discuss the results. We start with the LGW Hamiltonian in 
real space, 



f T 1 1 

H(m) = / d d f -g\Vm(f)\ 2 + -rm 2 (r) + um A (r) 

J \_Zi Zi 



(4.142) 



where u > 0. We first integrate out an infinitesimally thin shell of the thickness 
Sq in momentum space by taking 6 ~ 1. Then, 

Sq/A = 5A/A = 5\nA = b- 1 w lnfe (4.143) 

and we write 

m(f) = fh(r) + 5m(f) (4.144) 

where Sm(f) contains the Fourier components to be integrated out and rh the rest. 
After some approximations we get the following renormalization-group (differen- 
tial) equations for r and u [Huang, 1987]: 

dr(T) = 2r + 12A^fl 



dr V A 2 , 

= (4 - d)u - 36A d "V. (4.145) 



dr 

To analyse these equations, it is convenient to introduce the dimensionless 
coupling constants: 

r u 

x = a~^ y = A^ < 4 - 146 ) 



then, the RG equations become: 

drc 
d7 
dy 
dr 



2x + 12y(l -x) 

ey-36y 2 (4.147) 
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where 



e = 4 - d 



(4.148) 



We shall assume that e is small. The system has two fixed points at finite x and y, 
fixed by dx/dr = dy/dr = 0: 



The Gaussian fixed point is at x = y = whereas the "nontrivial" (also called 
the "Wilson-Fisher") fixed point lies in the upper (y > 0) half of the (x, y) plane 
only if e > 0, d < 4. In the lower half of the (x, y) plane, u is negative and the 
system is unstable. The nontrivial fixed point approaches the Gaussian fixed point 
as e — > 0. 

In the neighbourhood of the fixed points we linearize the RG equations 



x* = 
y*=0 

x* = -e/6 
V* = e/36 



Gaussian fixed point 



Nontrivial" fixed point. 



(4.149) 



x = x* + 5x, 



y = y* + Sy 



(4.150) 



and get the following linearized recursion relations: 



^Sx = (2 - I2y*)5x + 12(1 - x*)5y 
^-5y = (e-72y*)5y 



(4.151) 



which are written in the general form as: 




(4.152) 



At the Gaussian fixed point the linear matrix C is: 



2 12 
e 



(4.153) 



The matrix C has the following eigenvalues and eigenvectors: 
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A 2 = 6 ^(-(2-^/12) (4 - 154) 

The corresponding scaling fields are hi and h 2 . Under the (infinitesimal) RG 



transformations the scaling fields change as dhi/dr = \hi 

hi(r) = h i0 e x * T = h i0 b x > 



(4.155) 



Thus, the critical exponents are equal to the eigenvalues A^: x t = Ai = 2 and 

x 2 = A 2 = e 

In the neighbourhood of the nontrivial fixed point the linear transformation 
matrix C is: 



2(1 — 62/*) 12(1 -x* 
e - 72y* 

It has the following eigenvalues and eigenvectors: 



2(1 -e/6) 12(1 + e/6) 
-e 



(4.156) 



Ai 
A 2 



2 -e/3 



— e 



v 2 



6+e 



1+6/3 
1 



(4.157) 



Now we have the following picture. For d < 4 (e > 0) the nontrivial fixed 
point lies in the upper half of the (x, y) plane whereas the Gaussian fixed point 
is (as before) at x = y = 0, see Fig. 4.20. For the nontrivial fixed point, h 2 is 
irrelevant (it does not affect the critical exponents) and hi = r oc t the relevant 
scaling field whereas for the Gaussian fixed point, both fields are relevant. Thus, 
the line v 2 is the critical line. For d > 4 (e < 0)the nontrivial fixed point lies in the 
lower half of the (x, y) plane, i.e., in the unphysical region (u < 0), see Fig. 4.21, 
and the critical behaviour is governed by the Gaussian fixed point, d = 4 is the 
special case when both fixed points coincide. The linearized recursion relations 
give x 2 = 0, h 2 is marginal, and x 2 cannot tell us the direction of RG flow. In 
this case we must include higher-order terms which were neglected when we 
linearized the transformation matrix. These terms tell us that the scaling field h 2 
is irrelevant for u > at d = 4. 



4.4.3 Critical Exponents 

We have seen that the critical exponents of the nontrivial fixed point are 

x t = 2 -e/3 
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Figure 4.20: Fixed points and flow diagram of the Landau-Ginzburg-Wilson 
model for d < 4. The Gaussian fixed point at (x*,y*) = (0,0) has the mean- 
field exponents. The non-trivial fixed point at (x*,y*) = (— e/6, e/36) is unstable 
against the relevant perturbation 5x ~ t and stable against the irrelevant perturba- 
tion by. The system at A has t < and it flows towards the low-temperature fixed 
point, r becomes increasingly negative. At the point C the system has t > and 
it flows to the high-temperature fixed point, r becomes increasingly positive. The 
point B is on the critical line and the system flows towards the non-trivial fixed 
point. (After K. Huang, Statistical Mechanics) 
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yocu 




Figure 4.21: Fixed points and flow diagram of the Landau-Ginzburg- Wilson 
model for d > 4. The nontrivial fixed point is in the unphysical region. To 
stabilize the system in the lower half of the (x, y) plane, one needs an m 6 term. 
(After K. Huang, Statistical Mechanics) 



x h = l + d/2. (4.158) 

In fact, we calculated %h only for the Gaussian fixed point in the Eq. (4.138). 
However, the term with m 4 does not affect the value of Xh so that Xh of the non- 
trivial fixed point is the same. With these values of X{ we get from Eqs. (4.115, 
4.124, 4.125) to the first order in e for d < 4: 

e 
6 

1 _ e 

2 ~ 6 

1 + i 

3 + e 

1 e 

2 + 12 

0. (4.159) 



a = 

P = 

7 = 

5 = 

v = 

V = 
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For d > 4 (e < 0) we only have the Gaussian fixed point in the physical region 
(u > 0) with the exponents 

x t = 2 
x u e 

x h = l + d/2. (4.160) 
With Eqs. (4.115, 4.124, 4.125) we find: 

« d 
a = 2 "2 

(P = \~\) 
7 = 1 



(5 = 3 + e) 



v = 



1 

2 



r) = 0. (4.161) 

The above expressions for (3 and 5 are wrong, because we have neglected the 
m 4 term in considering the Gaussian fixed point. Below T c , r is negative and for 
m to be finite, the term urn 4 with u > must be included in the Hamiltonian. 
Although u is an irrelevant field, it helps to stabilize mat a finite value and in this 
way influences the critical behaviour. By treating the um 4 term in the Gaussian 
approximation (independent modes) we get the usual mean-field values for the 
exponents f3 and rj: 

5 = 3 (4.162) 

for all d > 4. These results are correct and also agree with the results of the 
Landau (MFA) model in Sec. 4.1.4. 

The Gaussian model gives the same critical exponents as the mean-field theory 
except for the specific-heat exponent a. Therefore the Gaussian fixed point is 
identified as the mean-field fixed point. In the Gaussian model the specific heat 
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diverges (a > 0) for d < 4. This divergence is caused by the long-wavelength 
fluctuations which were neglected in the calculation of the free energy in the MF 
approximation. 

In d = 3 (e = 1) we find a xs 0.17 and v m 0.58. The results of this paragraph 
are based on an expansion of d around d — 4 and are valid for small e > or for 
d > 4. Nevertheless, we put e = 1 when d — 3. That means that e is not small and 
higher-order terms in e must be included. This is then the e— expansion, which 
gives good values for the critical exponents. 
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Appendix A 
Distribution Functions 



We have to distinguish three kinds of systems: 

(a) identical but distinguishable particles ("classical particles") 

(b) identical indistinguishable particles with half-integer spin (Fermi-ons) 

(c) identical indistinguishable particles with integer spins (Bosons) 
The occupation probability of a state % with energy E { is: 

(a-Maxwell-Boltzmann): 




(A.1) 



2 



Maxwell Boltzmann distrubution 



T = 1 000 K 



g 



2: 



T = 5 000 K 












0.5 



1.5 



2 



Energy (eV) 



(b-Fermi-Dirac): 




(A.2) 



(c-Bose-Einstein): 




(A.3) 
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Fermi Dirac distrubution 




= 1 eV 
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Appendix B 
Maxwell Relations 



Let df be a total differential of the function f(x,y). Then, 

df = u(x, y)dx + v (x, y)dy (B.l) 

with 



and 

'"-a <*r\ = (m = (¥\. (b.3) 



\dy I \ / \ &r<9y / \dx I 

This is a very useful relation that - when applied to thermodynamic functions - 
gives the Maxwell relations: 



'&T\ _(dp\ 
fdT\ fdV\ 

[dp) s -[ds); 

'dS\ 



(B.4) 

(B.5) 
(B.6) 
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Appendix C 



Basic Thermodynamic Relations of 
Magnetic Systems Revisited 

The internal energy is: 

u(s,M) = -Y t JiASiS i ) (c.i) 

(ij) 

where S is the entropy, M is the total magnetic moment (in dimensionless units), 
M = Nm, m is the magnetic moment per lattice site, m = (Si). J it j is the 
exchange interaction between the spins % and j. The (magnetic) enthalpy is (here 
we will use E for the enthalpy, to avoid confusion with the field H): 

E(S, H) — U — EM = (H) (C.2) 

where H is the magnetic field in units of energy. The Helmholtz free energy is: 

F(T, M) = —k B T In Z(T, M, N) = U — TS (C.3) 

and its total differential is 

dF = -SdT + HdM (C.4) 

(see Table II of Part I; the number of spins N is constant): The Gibbs free energy 
is 

G(T, H) = —k B T In Y(T, H) — E — TS 

= F - EM = U — TS — EM (C.5) 



135 



I. Vilfan 



Statistical Mechanics 



Appendix C 



and its differential is 



dG = -SdT-MdH 



(C6) 



Here, we introduced also the corresponding partition functions Z and Y. From 
Eq. (C.4) we get immediately the equation of state: 



Of course, the equations (C.7) and (C.8) must lead to the same equation, they are 
just two ways of writing the equation of state. When we want to study systems at 
constant H (and not at constant M ), the more convenient thermodynamic potential 
is the Gibbsfree energy G(T, H, N) instead of the Helmholtz free energy. 
Other useful thermodynamic relations are: 




(C.7) 



From Eq. (C.6) we also get an equation of state: 




(C.8) 




(C9) 




(C.IO) 




(Cll) 




(C12) 
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